1 #ifdef RCSID
2 static char RCSid[] =
3 "$Header$";
4 #endif
5
6 /*
7 * Copyright (c) 2000, 2002 Michael J. Roberts. All Rights Reserved.
8 *
9 * Please see the accompanying license file, LICENSE.TXT, for information
10 * on using and copying this software.
11 */
12 /*
13 Name
14 vmbignum.cpp - big number metaclass
15 Function
16
17 Notes
18
19 Modified
20 02/18/00 MJRoberts - Creation
21 */
22
23 #include <stdlib.h>
24 #include <string.h>
25 #include <assert.h>
26 #include <limits.h>
27 #include <stdarg.h>
28
29 #include "t3std.h"
30 #include "vmtype.h"
31 #include "vmmcreg.h"
32 #include "vmbignum.h"
33 #include "vmobj.h"
34 #include "vmerr.h"
35 #include "vmerrnum.h"
36 #include "vmfile.h"
37 #include "vmpool.h"
38 #include "vmstack.h"
39 #include "utf8.h"
40 #include "vmstr.h"
41 #include "vmbif.h"
42 #include "vmmeta.h"
43 #include "vmlst.h"
44
45
46 /* ------------------------------------------------------------------------ */
47 /*
48 * statics
49 */
50
51 /* metaclass registration object */
52 static CVmMetaclassBigNum metaclass_reg_obj;
53 CVmMetaclass *CVmObjBigNum::metaclass_reg_ = &metaclass_reg_obj;
54
55 /* function table */
56 int (CVmObjBigNum::
57 *CVmObjBigNum::func_table_[])(VMG_ vm_obj_id_t self,
58 vm_val_t *retval, uint *argc) =
59 {
60 &CVmObjBigNum::getp_undef,
61 &CVmObjBigNum::getp_format,
62 &CVmObjBigNum::getp_equal_rnd,
63 &CVmObjBigNum::getp_get_prec,
64 &CVmObjBigNum::getp_set_prec,
65 &CVmObjBigNum::getp_frac,
66 &CVmObjBigNum::getp_whole,
67 &CVmObjBigNum::getp_round_dec,
68 &CVmObjBigNum::getp_abs,
69 &CVmObjBigNum::getp_ceil,
70 &CVmObjBigNum::getp_floor,
71 &CVmObjBigNum::getp_get_scale,
72 &CVmObjBigNum::getp_scale,
73 &CVmObjBigNum::getp_negate,
74 &CVmObjBigNum::getp_copy_sign,
75 &CVmObjBigNum::getp_is_neg,
76 &CVmObjBigNum::getp_remainder,
77 &CVmObjBigNum::getp_sin,
78 &CVmObjBigNum::getp_cos,
79 &CVmObjBigNum::getp_tan,
80 &CVmObjBigNum::getp_deg2rad,
81 &CVmObjBigNum::getp_rad2deg,
82 &CVmObjBigNum::getp_asin,
83 &CVmObjBigNum::getp_acos,
84 &CVmObjBigNum::getp_atan,
85 &CVmObjBigNum::getp_sqrt,
86 &CVmObjBigNum::getp_ln,
87 &CVmObjBigNum::getp_exp,
88 &CVmObjBigNum::getp_log10,
89 &CVmObjBigNum::getp_pow,
90 &CVmObjBigNum::getp_sinh,
91 &CVmObjBigNum::getp_cosh,
92 &CVmObjBigNum::getp_tanh,
93 &CVmObjBigNum::getp_pi,
94 &CVmObjBigNum::getp_e
95 };
96
97
98 /* constant value 1 */
99 const unsigned char CVmObjBigNum::one_[] =
100 {
101 /* number of digits - we just need one to represent the integer 1 */
102 0x01, 0x00,
103
104 /* exponent */
105 0x01, 0x00,
106
107 /* flags */
108 0x00,
109
110 /* mantissa - just one digit is needed, but pad out the byte with zero */
111 0x10
112 };
113
114 #if 0
115 /*
116 * Pi to 2048 digits
117 */
118 const unsigned char CVmObjBigNum::pi_[] =
119 {
120 /* number of digits of precision */
121 0x00, 0x08,
122
123 /* base-10 exponent */
124 0x01, 0x00,
125
126 /* flags */
127 0x00,
128
129 /*
130 * the first 2048 decimal digits of pi, packed two to a byte (typed
131 * in from memory, I hope I got everything right :)
132 */
133
134 /* 1-100 */
135 0x31, 0x41, 0x59, 0x26, 0x53, 0x58, 0x97, 0x93, 0x23, 0x84,
136 0x62, 0x64, 0x33, 0x83, 0x27, 0x95, 0x02, 0x88, 0x41, 0x97,
137 0x16, 0x93, 0x99, 0x37, 0x51, 0x05, 0x82, 0x09, 0x74, 0x94,
138 0x45, 0x92, 0x30, 0x78, 0x16, 0x40, 0x62, 0x86, 0x20, 0x89,
139 0x98, 0x62, 0x80, 0x34, 0x82, 0x53, 0x42, 0x11, 0x70, 0x67,
140
141 /* 101-200 */
142 0x98, 0x21, 0x48, 0x08, 0x65, 0x13, 0x28, 0x23, 0x06, 0x64,
143 0x70, 0x93, 0x84, 0x46, 0x09, 0x55, 0x05, 0x82, 0x23, 0x17,
144 0x25, 0x35, 0x94, 0x08, 0x12, 0x84, 0x81, 0x11, 0x74, 0x50,
145 0x28, 0x41, 0x02, 0x70, 0x19, 0x38, 0x52, 0x11, 0x05, 0x55,
146 0x96, 0x44, 0x62, 0x29, 0x48, 0x95, 0x49, 0x30, 0x38, 0x19,
147
148 /* 201-300 */
149 0x64, 0x42, 0x88, 0x10, 0x97, 0x56, 0x65, 0x93, 0x34, 0x46,
150 0x12, 0x84, 0x75, 0x64, 0x82, 0x33, 0x78, 0x67, 0x83, 0x16,
151 0x52, 0x71, 0x20, 0x19, 0x09, 0x14, 0x56, 0x48, 0x56, 0x69,
152 0x23, 0x46, 0x03, 0x48, 0x61, 0x04, 0x54, 0x32, 0x66, 0x48,
153 0x21, 0x33, 0x93, 0x60, 0x72, 0x60, 0x24, 0x91, 0x41, 0x27,
154
155 /* 301-400 */
156 0x37, 0x24, 0x58, 0x70, 0x06, 0x60, 0x63, 0x15, 0x58, 0x81,
157 0x74, 0x88, 0x15, 0x20, 0x92, 0x09, 0x62, 0x82, 0x92, 0x54,
158 0x09, 0x17, 0x15, 0x36, 0x43, 0x67, 0x89, 0x25, 0x90, 0x36,
159 0x00, 0x11, 0x33, 0x05, 0x30, 0x54, 0x88, 0x20, 0x46, 0x65,
160 0x21, 0x38, 0x41, 0x46, 0x95, 0x19, 0x41, 0x51, 0x16, 0x09,
161
162 /* 401-500 */
163 0x43, 0x30, 0x57, 0x27, 0x03, 0x65, 0x75, 0x95, 0x91, 0x95,
164 0x30, 0x92, 0x18, 0x61, 0x17, 0x38, 0x19, 0x32, 0x61, 0x17,
165 0x93, 0x10, 0x51, 0x18, 0x54, 0x80, 0x74, 0x46, 0x23, 0x79,
166 0x96, 0x27, 0x49, 0x56, 0x73, 0x51, 0x88, 0x57, 0x52, 0x72,
167 0x48, 0x91, 0x22, 0x79, 0x38, 0x18, 0x30, 0x11, 0x94, 0x91,
168
169 /* 501-600 */
170 0x29, 0x83, 0x36, 0x73, 0x36, 0x24, 0x40, 0x65, 0x66, 0x43,
171 0x08, 0x60, 0x21, 0x39, 0x49, 0x46, 0x39, 0x52, 0x24, 0x73,
172 0x71, 0x90, 0x70, 0x21, 0x79, 0x86, 0x09, 0x43, 0x70, 0x27,
173 0x70, 0x53, 0x92, 0x17, 0x17, 0x62, 0x93, 0x17, 0x67, 0x52,
174 0x38, 0x46, 0x74, 0x81, 0x84, 0x67, 0x66, 0x94, 0x05, 0x13,
175
176 /* 601-700 */
177 0x20, 0x00, 0x56, 0x81, 0x27, 0x14, 0x52, 0x63, 0x56, 0x08,
178 0x27, 0x78, 0x57, 0x71, 0x34, 0x27, 0x57, 0x78, 0x96, 0x09,
179 0x17, 0x36, 0x37, 0x17, 0x87, 0x21, 0x46, 0x84, 0x40, 0x90,
180 0x12, 0x24, 0x95, 0x34, 0x30, 0x14, 0x65, 0x49, 0x58, 0x53,
181 0x71, 0x05, 0x07, 0x92, 0x27, 0x96, 0x89, 0x25, 0x89, 0x23,
182
183 /* 701-800 */
184 0x54, 0x20, 0x19, 0x95, 0x61, 0x12, 0x12, 0x90, 0x21, 0x96,
185 0x08, 0x64, 0x03, 0x44, 0x18, 0x15, 0x98, 0x13, 0x62, 0x97,
186 0x74, 0x77, 0x13, 0x09, 0x96, 0x05, 0x18, 0x70, 0x72, 0x11,
187 0x34, 0x99, 0x99, 0x99, 0x83, 0x72, 0x97, 0x80, 0x49, 0x95,
188 0x10, 0x59, 0x73, 0x17, 0x32, 0x81, 0x60, 0x96, 0x31, 0x85,
189
190 /* 801-900 */
191 0x95, 0x02, 0x44, 0x59, 0x45, 0x53, 0x46, 0x90, 0x83, 0x02,
192 0x64, 0x25, 0x22, 0x30, 0x82, 0x53, 0x34, 0x46, 0x85, 0x03,
193 0x52, 0x61, 0x93, 0x11, 0x88, 0x17, 0x10, 0x10, 0x00, 0x31,
194 0x37, 0x83, 0x87, 0x52, 0x88, 0x65, 0x87, 0x53, 0x32, 0x08,
195 0x38, 0x14, 0x20, 0x61, 0x71, 0x77, 0x66, 0x91, 0x47, 0x30,
196
197 /* 901-1000 */
198 0x35, 0x98, 0x25, 0x34, 0x90, 0x42, 0x87, 0x55, 0x46, 0x87,
199 0x31, 0x15, 0x95, 0x62, 0x86, 0x38, 0x82, 0x35, 0x37, 0x87,
200 0x59, 0x37, 0x51, 0x95, 0x77, 0x81, 0x85, 0x77, 0x80, 0x53,
201 0x21, 0x71, 0x22, 0x68, 0x06, 0x61, 0x30, 0x01, 0x92, 0x78,
202 0x76, 0x61, 0x11, 0x95, 0x90, 0x92, 0x16, 0x42, 0x01, 0x98,
203
204 /* 1001-1100 */
205 0x93, 0x80, 0x95, 0x25, 0x72, 0x01, 0x06, 0x54, 0x85, 0x86,
206 0x32, 0x78, 0x86, 0x59, 0x36, 0x15, 0x33, 0x81, 0x82, 0x79,
207 0x68, 0x23, 0x03, 0x01, 0x95, 0x20, 0x35, 0x30, 0x18, 0x52,
208 0x96, 0x89, 0x95, 0x77, 0x36, 0x22, 0x59, 0x94, 0x13, 0x89,
209 0x12, 0x49, 0x72, 0x17, 0x75, 0x28, 0x34, 0x79, 0x13, 0x15,
210
211 /* 1101-1200 */
212 0x15, 0x57, 0x48, 0x57, 0x24, 0x24, 0x54, 0x15, 0x06, 0x95,
213 0x95, 0x08, 0x29, 0x53, 0x31, 0x16, 0x86, 0x17, 0x27, 0x85,
214 0x58, 0x89, 0x07, 0x50, 0x98, 0x38, 0x17, 0x54, 0x63, 0x74,
215 0x64, 0x93, 0x93, 0x19, 0x25, 0x50, 0x60, 0x40, 0x09, 0x27,
216 0x70, 0x16, 0x71, 0x13, 0x90, 0x09, 0x84, 0x88, 0x24, 0x01,
217
218 /* 1201-1300 */
219 0x28, 0x58, 0x36, 0x16, 0x03, 0x56, 0x37, 0x07, 0x66, 0x01,
220 0x04, 0x71, 0x01, 0x81, 0x94, 0x29, 0x55, 0x59, 0x61, 0x98,
221 0x94, 0x67, 0x67, 0x83, 0x74, 0x49, 0x44, 0x82, 0x55, 0x37,
222 0x97, 0x74, 0x72, 0x68, 0x47, 0x10, 0x40, 0x47, 0x53, 0x46,
223 0x46, 0x20, 0x80, 0x46, 0x68, 0x42, 0x59, 0x06, 0x94, 0x91,
224
225 /* 1301-1400 */
226 0x29, 0x33, 0x13, 0x67, 0x70, 0x28, 0x98, 0x91, 0x52, 0x10,
227 0x47, 0x52, 0x16, 0x20, 0x56, 0x96, 0x60, 0x24, 0x05, 0x80,
228 0x38, 0x15, 0x01, 0x93, 0x51, 0x12, 0x53, 0x38, 0x24, 0x30,
229 0x03, 0x55, 0x87, 0x64, 0x02, 0x47, 0x49, 0x64, 0x73, 0x26,
230 0x39, 0x14, 0x19, 0x92, 0x72, 0x60, 0x42, 0x69, 0x92, 0x27,
231
232 /* 1401-1500 */
233 0x96, 0x78, 0x23, 0x54, 0x78, 0x16, 0x36, 0x00, 0x93, 0x41,
234 0x72, 0x16, 0x41, 0x21, 0x99, 0x24, 0x58, 0x63, 0x15, 0x03,
235 0x02, 0x86, 0x18, 0x29, 0x74, 0x55, 0x57, 0x06, 0x74, 0x98,
236 0x38, 0x50, 0x54, 0x94, 0x58, 0x85, 0x86, 0x92, 0x69, 0x95,
237 0x69, 0x09, 0x27, 0x21, 0x07, 0x97, 0x50, 0x93, 0x02, 0x95,
238
239 /* 1501-1600 */
240 0x53, 0x21, 0x16, 0x53, 0x44, 0x98, 0x72, 0x02, 0x75, 0x59,
241 0x60, 0x23, 0x64, 0x80, 0x66, 0x54, 0x99, 0x11, 0x98, 0x81,
242 0x83, 0x47, 0x97, 0x75, 0x35, 0x66, 0x36, 0x98, 0x07, 0x42,
243 0x65, 0x42, 0x52, 0x78, 0x62, 0x55, 0x18, 0x18, 0x41, 0x75,
244 0x74, 0x67, 0x28, 0x90, 0x97, 0x77, 0x72, 0x79, 0x38, 0x00,
245
246 /* 1601-1700 */
247 0x08, 0x16, 0x47, 0x06, 0x00, 0x16, 0x14, 0x52, 0x49, 0x19,
248 0x21, 0x73, 0x21, 0x72, 0x14, 0x77, 0x23, 0x50, 0x14, 0x14,
249 0x41, 0x97, 0x35, 0x68, 0x54, 0x81, 0x61, 0x36, 0x11, 0x57,
250 0x35, 0x25, 0x52, 0x13, 0x34, 0x75, 0x74, 0x18, 0x49, 0x46,
251 0x84, 0x38, 0x52, 0x33, 0x23, 0x90, 0x73, 0x94, 0x14, 0x33,
252
253 /* 1701-1800 */
254 0x34, 0x54, 0x77, 0x62, 0x41, 0x68, 0x62, 0x51, 0x89, 0x83,
255 0x56, 0x94, 0x85, 0x56, 0x20, 0x99, 0x21, 0x92, 0x22, 0x18,
256 0x42, 0x72, 0x55, 0x02, 0x54, 0x25, 0x68, 0x87, 0x67, 0x17,
257 0x90, 0x49, 0x46, 0x01, 0x65, 0x34, 0x66, 0x80, 0x49, 0x88,
258 0x62, 0x72, 0x32, 0x79, 0x17, 0x86, 0x08, 0x57, 0x84, 0x38,
259
260 /* 1801-1900 */
261 0x38, 0x27, 0x96, 0x79, 0x76, 0x68, 0x14, 0x54, 0x10, 0x09,
262 0x53, 0x88, 0x37, 0x86, 0x36, 0x09, 0x50, 0x68, 0x00, 0x64,
263 0x22, 0x51, 0x25, 0x20, 0x51, 0x17, 0x39, 0x29, 0x84, 0x89,
264 0x60, 0x84, 0x12, 0x84, 0x88, 0x62, 0x69, 0x45, 0x60, 0x42,
265 0x41, 0x96, 0x52, 0x85, 0x02, 0x22, 0x10, 0x66, 0x11, 0x86,
266
267 /* 1901-2000 */
268 0x30, 0x67, 0x44, 0x27, 0x86, 0x22, 0x03, 0x91, 0x94, 0x94,
269 0x50, 0x47, 0x12, 0x37, 0x13, 0x78, 0x69, 0x60, 0x95, 0x63,
270 0x64, 0x37, 0x19, 0x17, 0x28, 0x74, 0x67, 0x76, 0x46, 0x57,
271 0x57, 0x39, 0x62, 0x41, 0x38, 0x90, 0x86, 0x58, 0x32, 0x64,
272 0x59, 0x95, 0x81, 0x33, 0x90, 0x47, 0x80, 0x27, 0x59, 0x00,
273
274 /* 2001-2048 */
275 0x99, 0x46, 0x57, 0x64, 0x07, 0x89, 0x51, 0x26, 0x94, 0x68,
276 0x39, 0x83, 0x52, 0x59, 0x57, 0x09, 0x82, 0x58, 0x22, 0x62,
277 0x05, 0x22, 0x48, 0x94
278 };
279 #endif
280
281 /* ------------------------------------------------------------------------ */
282 /*
283 * Static creation methods. These routines allocate an object ID and
284 * create a new object.
285 */
286
287
288 /* create dynamically using stack arguments */
create_from_stack(VMG_ const uchar ** pc_ptr,uint argc)289 vm_obj_id_t CVmObjBigNum::create_from_stack(VMG_ const uchar **pc_ptr,
290 uint argc)
291 {
292 vm_obj_id_t id;
293 vm_val_t *val;
294 size_t digits;
295 const char *strval = 0;
296 const CVmObjBigNum *objval = 0;
297
298 /* check arguments */
299 if (argc < 1 || argc > 2)
300 err_throw(VMERR_WRONG_NUM_OF_ARGS);
301
302 /*
303 * check the first argument, which gives the source value - this can be
304 * an integer, a string, or another BigNumber value
305 */
306 val = G_stk->get(0);
307 if (val->typ == VM_INT)
308 {
309 /* we'll use the integer value */
310 }
311 else if (val->typ == VM_OBJ && is_bignum_obj(vmg_ val->val.obj))
312 {
313 /* we'll use the BigNumber value as the source */
314 objval = (CVmObjBigNum *)vm_objp(vmg_ val->val.obj);
315 }
316 else if ((strval = val->get_as_string(vmg0_)) != 0)
317 {
318 /* we'll parse the string value */
319 }
320 else
321 {
322 /* it's not a source type we accept */
323 err_throw(VMERR_NUM_VAL_REQD);
324 }
325
326 /*
327 * get the precision value, if provided; if not, infer it from the
328 * source value
329 */
330 if (argc > 1)
331 {
332 /* make sure it's an integer */
333 if (G_stk->get(1)->typ != VM_INT)
334 err_throw(VMERR_INT_VAL_REQD);
335
336 /* retrieve the value */
337 digits = (size_t)G_stk->get(1)->val.intval;
338 }
339 else if (strval != 0)
340 {
341 utf8_ptr p;
342 size_t rem;
343 size_t prec;
344 int pt;
345 int significant;
346 int digcnt;
347
348 /* set up to scan the string */
349 p.set((char *)strval + VMB_LEN);
350 rem = vmb_get_len(strval);
351
352 /* skip leading spaces */
353 for ( ; rem != 0 && is_space(p.getch()) ; p.inc(&rem)) ;
354
355 /* skip the sign if present */
356 if (rem != 0 && (p.getch() == '-' || p.getch() == '+'))
357 p.inc(&rem);
358
359 /* we haven't yet found a significant leading digit */
360 significant = FALSE;
361
362 /* scan digits */
363 for (prec = 0, digcnt = 0, pt = FALSE ; rem != 0 ; p.inc(&rem))
364 {
365 wchar_t ch;
366
367 /* get this character */
368 ch = p.getch();
369
370 /* see what we have */
371 if (is_digit(ch))
372 {
373 /*
374 * if it's not a zero, note that we've found a
375 * significant digit
376 */
377 if (ch != '0')
378 significant = TRUE;
379
380 /*
381 * if we have found a significant digit so far, count
382 * this one - do not count leading zeroes, whether they
383 * occur before or after the decimal point
384 */
385 if (significant)
386 ++prec;
387
388 /* count the digit in any case */
389 ++digcnt;
390 }
391 else if (ch == '.' && !pt)
392 {
393 /* decimal point - note it and keep going */
394 pt = TRUE;
395 }
396 else if (ch == 'e' || ch == 'E')
397 {
398 /* exponent - that's the end of the mantissa */
399 break;
400 }
401 else
402 {
403 /* we've reached the end of the number */
404 break;
405 }
406 }
407
408 /*
409 * if the precision is zero, the number must be zero - use the
410 * actual number of digits for the default precision, so that a
411 * value specified as "0.0000000" has eight digits of precision
412 */
413 if (prec == 0)
414 prec = digcnt;
415
416 /* use the precision necessary to store the entire string */
417 digits = prec;
418 }
419 else if (objval != 0)
420 {
421 /* use the same precision as the source BigNumber value */
422 digits = get_prec(objval->get_ext());
423 }
424 else
425 {
426 /* use a default precision */
427 digits = 32;
428 }
429
430 /* create the value */
431 if (strval != 0)
432 {
433 /* create the value from the string */
434 id = vm_new_id(vmg_ FALSE, FALSE, FALSE);
435 new (vmg_ id) CVmObjBigNum(vmg_ strval + VMB_LEN,
436 vmb_get_len(strval), digits);
437 }
438 else if (objval != 0)
439 {
440 vm_val_t new_val;
441
442 /* create the value based on the BigNumber value */
443 round_val(vmg_ &new_val, objval->get_ext(), digits, TRUE);
444
445 /* return the new object ID */
446 id = new_val.val.obj;
447 }
448 else
449 {
450 /* create the value based on the integer value */
451 id = vm_new_id(vmg_ FALSE, FALSE, FALSE);
452 new (vmg_ id) CVmObjBigNum(vmg_ val->val.intval, digits);
453 }
454
455 /* discard arguments */
456 G_stk->discard(argc);
457
458 /* return the new object */
459 return id;
460 }
461
462 /* create with a given precision */
create(VMG_ int in_root_set,size_t digits)463 vm_obj_id_t CVmObjBigNum::create(VMG_ int in_root_set, size_t digits)
464 {
465 vm_obj_id_t id = vm_new_id(vmg_ in_root_set, FALSE, FALSE);
466 new (vmg_ id) CVmObjBigNum(vmg_ digits);
467 return id;
468 }
469
470 /* create from an integer value */
create(VMG_ int in_root_set,long val,size_t digits)471 vm_obj_id_t CVmObjBigNum::create(VMG_ int in_root_set,
472 long val, size_t digits)
473 {
474 vm_obj_id_t id = vm_new_id(vmg_ in_root_set, FALSE, FALSE);
475 new (vmg_ id) CVmObjBigNum(vmg_ val, digits);
476 return id;
477 }
478
479 /* ------------------------------------------------------------------------ */
480 /*
481 * Constructors. These are called indirectly through our static
482 * creation methods.
483 */
484
485 /*
486 * Create with no extension
487 */
CVmObjBigNum()488 CVmObjBigNum::CVmObjBigNum()
489 {
490 /* no extension */
491 ext_ = 0;
492 }
493
494 /*
495 * Create with a given precision
496 */
CVmObjBigNum(VMG_ size_t digits)497 CVmObjBigNum::CVmObjBigNum(VMG_ size_t digits)
498 {
499 /* allocate space */
500 alloc_bignum(vmg_ digits);
501 }
502
503 /*
504 * Create with a given integer value
505 */
CVmObjBigNum(VMG_ long val,size_t digits)506 CVmObjBigNum::CVmObjBigNum(VMG_ long val, size_t digits)
507 {
508 /* allocate space */
509 alloc_bignum(vmg_ digits);
510
511 /* set the value */
512 set_int_val(val);
513 }
514
515 /*
516 * Create with a given string as the source value
517 */
CVmObjBigNum(VMG_ const char * str,size_t len,size_t digits)518 CVmObjBigNum::CVmObjBigNum(VMG_ const char *str, size_t len, size_t digits)
519 {
520 /* allocate space */
521 alloc_bignum(vmg_ digits);
522
523 /* set the value */
524 set_str_val(str, len);
525 }
526
527 /* ------------------------------------------------------------------------ */
528 /*
529 * Delete
530 */
notify_delete(VMG_ int in_root_set)531 void CVmObjBigNum::notify_delete(VMG_ int in_root_set)
532 {
533 /*
534 * free our extension - do this only if it's not in the root set,
535 * because extension will be directly in the image data for a root
536 * set object
537 */
538 if (ext_ != 0 && !in_root_set)
539 G_mem->get_var_heap()->free_mem(ext_);
540 }
541
542
543 /* ------------------------------------------------------------------------ */
544 /*
545 * Allocate space for a given precision
546 */
alloc_bignum(VMG_ size_t digits)547 void CVmObjBigNum::alloc_bignum(VMG_ size_t digits)
548 {
549 /* allocate space for the given number of elements */
550 ext_ = (char *)G_mem->get_var_heap()
551 ->alloc_mem(calc_alloc(digits), this);
552
553 /* set the precision */
554 set_prec(ext_, digits);
555
556 /* initialize the value to zero */
557 set_int_val(0);
558
559 /* clear the flags */
560 ext_[VMBN_FLAGS] = 0;
561 }
562
563 /*
564 * Calculate the amount of memory we need for a given number of digits
565 * of precision
566 */
calc_alloc(size_t digits)567 size_t CVmObjBigNum::calc_alloc(size_t digits)
568 {
569 /*
570 * we need the header (UINT2, INT2, BYTE), plus one byte for each
571 * two decimal digits
572 */
573 return (2 + 2 + 1) + ((digits + 1)/2);
574 }
575
576
577 /* ------------------------------------------------------------------------ */
578 /*
579 * Write to a 'data' mode file
580 */
write_to_data_file(osfildef * fp)581 int CVmObjBigNum::write_to_data_file(osfildef *fp)
582 {
583 char buf[16];
584
585 /* write the number of digits (i.e., the precision) */
586 oswp2(buf, get_prec(ext_));
587 if (osfwb(fp, buf, 2))
588 return 1;
589
590 /* write our entire extension */
591 if (osfwb(fp, ext_, calc_alloc(get_prec(ext_))))
592 return 1;
593
594 /* success */
595 return 0;
596 }
597
598 /*
599 * Read from a 'data' mode file and instantiate a new BigNumber object to
600 * hold the result
601 */
read_from_data_file(VMG_ vm_val_t * retval,osfildef * fp)602 int CVmObjBigNum::read_from_data_file(VMG_ vm_val_t *retval, osfildef *fp)
603 {
604 char buf[16];
605 size_t prec;
606 CVmObjBigNum *bignum;
607
608 /* read the precision */
609 if (osfrb(fp, buf, 2))
610 return 1;
611 prec = osrp2(buf);
612
613 /* create a BigNumber with the required precision */
614 retval->set_obj(create(vmg_ FALSE, prec));
615 bignum = (CVmObjBigNum *)vm_objp(vmg_ retval->val.obj);
616
617 /* read the bytes into the new object's extension */
618 if (osfrb(fp, bignum->get_ext(), calc_alloc(prec)))
619 return 1;
620
621 /* success */
622 return 0;
623 }
624
625
626 /* ------------------------------------------------------------------------ */
627 /*
628 * Set my value to a given integer value
629 */
set_int_val(long val)630 void CVmObjBigNum::set_int_val(long val)
631 {
632 size_t prec;
633 int exp;
634
635 /* get the precision */
636 prec = get_prec(ext_);
637
638 /* set the type to number */
639 set_type(ext_, VMBN_T_NUM);
640
641 /* set the sign bit */
642 if (val < 0)
643 {
644 /* set the value negative */
645 set_neg(ext_, TRUE);
646
647 /* use the absolute value for the mantissa */
648 val = -val;
649 }
650 else
651 {
652 /* set the value positive */
653 set_neg(ext_, FALSE);
654 }
655
656 /* initially zero the mantissa */
657 memset(ext_ + VMBN_MANT, 0, (prec + 1)/2);
658
659 /* initialize the exponent to zero */
660 exp = 0;
661
662 /* shift the integer value into the value */
663 while (val != 0)
664 {
665 unsigned int dig;
666
667 /* get the low-order digit of the value */
668 dig = (unsigned int)(val % 10);
669
670 /* shift the value one place */
671 val /= 10;
672
673 /*
674 * shift our number one place right to accommodate this next
675 * higher-order digit
676 */
677 shift_right(ext_, 1);
678
679 /* set the new most significant digit */
680 set_dig(ext_, 0, dig);
681
682 /* that's another factor of 10 */
683 ++exp;
684 }
685
686 /* set the exponent */
687 set_exp(ext_, exp);
688
689 /* normalize the number */
690 normalize(ext_);
691 }
692
693 /*
694 * Set my value to a string
695 */
set_str_val(const char * str,size_t len)696 void CVmObjBigNum::set_str_val(const char *str, size_t len)
697 {
698 size_t prec;
699 int exp;
700 utf8_ptr p;
701 size_t rem;
702 int neg;
703 size_t idx;
704 int pt;
705 int significant;
706
707 /* get the precision */
708 prec = get_prec(ext_);
709
710 /* set the type to number */
711 set_type(ext_, VMBN_T_NUM);
712
713 /* initially zero the mantissa */
714 memset(ext_ + VMBN_MANT, 0, (prec + 1)/2);
715
716 /* set up to scan the string */
717 p.set((char *)str);
718 rem = len;
719
720 /* initialize the exponent to zero */
721 exp = 0;
722
723 /* presume it will be positive */
724 neg = FALSE;
725
726 /* skip leading spaces */
727 while (rem != 0 && is_space(p.getch()))
728 p.inc(&rem);
729
730 /* check for a sign */
731 if (rem != 0 && p.getch() == '+')
732 {
733 /* skip the sign */
734 p.inc(&rem);
735 }
736 else if (rem != 0 && p.getch() == '-')
737 {
738 /* note the sign and skip it */
739 neg = TRUE;
740 p.inc(&rem);
741 }
742
743 /* set the sign */
744 set_neg(ext_, neg);
745
746 /* we haven't yet found a significant digit */
747 significant = FALSE;
748
749 /* shift the digits into the value */
750 for (idx = 0, pt = FALSE ; rem != 0 ; p.inc(&rem))
751 {
752 wchar_t ch;
753
754 /* get this character */
755 ch = p.getch();
756
757 /* check for a digit */
758 if (is_digit(ch))
759 {
760 /*
761 * if it's not a zero, we're definitely into the significant
762 * part of the number
763 */
764 if (ch != '0')
765 significant = TRUE;
766
767 /*
768 * if it's significant, add it to the number - skip leading
769 * zeroes, since they add no information to the number
770 */
771 if (significant)
772 {
773 /* if we're not out of precision, add the digit */
774 if (idx < prec)
775 {
776 /* set the next digit */
777 set_dig(ext_, idx, value_of_digit(ch));
778
779 /* move on to the next digit position */
780 ++idx;
781 }
782
783 /*
784 * that's another factor of 10 if we haven't found the
785 * decimal point (whether or not we're out of precision
786 * to actually store the digit, count the increase in
787 * the exponent)
788 */
789 if (!pt)
790 ++exp;
791 }
792 else if (pt)
793 {
794 /*
795 * we haven't yet found a significant digit, so this is
796 * a leading zero, but we have found the decimal point -
797 * this reduces the exponent by one
798 */
799 --exp;
800 }
801 }
802 else if (ch == '.' && !pt)
803 {
804 /*
805 * this is the decimal point - note it; from now on, we
806 * won't increase the exponent as we add digits
807 */
808 pt = TRUE;
809 }
810 else if (ch == 'e' || ch == 'E')
811 {
812 int acc;
813 int exp_neg = FALSE;
814
815 /* exponent - skip the 'e' */
816 p.inc(&rem);
817
818 /* check for a sign */
819 if (rem != 0 && p.getch() == '+')
820 {
821 /* skip the sign */
822 p.inc(&rem);
823 }
824 else if (rem != 0 && p.getch() == '-')
825 {
826 /* skip it and note it */
827 p.inc(&rem);
828 exp_neg = TRUE;
829 }
830
831 /* parse the exponent */
832 for (acc = 0 ; rem != 0 ; p.inc(&rem))
833 {
834 wchar_t ch;
835
836 /* if this is a digit, add it to the exponent */
837 ch = p.getch();
838 if (is_digit(ch))
839 {
840 /* accumulate the digit */
841 acc *= 10;
842 acc += value_of_digit(ch);
843 }
844 else
845 {
846 /* that's it for the exponent */
847 break;
848 }
849 }
850
851 /* if it's negative, apply the sign */
852 if (exp_neg)
853 acc = -acc;
854
855 /* add the exponent to the one we've calculated */
856 exp += acc;
857
858 /* after the exponent, we're done with the whole number */
859 break;
860 }
861 else
862 {
863 /* it looks like we've reached the end of the number */
864 break;
865 }
866 }
867
868 /* set the exponent */
869 set_exp(ext_, exp);
870
871 /* normalize the number */
872 normalize(ext_);
873 }
874
875
876 /* ------------------------------------------------------------------------ */
877 /*
878 * Convert to an integer value
879 */
convert_to_int()880 long CVmObjBigNum::convert_to_int()
881 {
882 const long long_max_over_10 = LONG_MAX/10;
883 const long long_min_over_10 = LONG_MIN/10;
884 size_t prec = get_prec(ext_);
885 int is_neg = get_neg(ext_);
886 int exp = get_exp(ext_);
887 size_t idx;
888 size_t stop_idx;
889 long acc;
890 int round_inc;
891
892 /* start the accumulator at zero */
893 acc = 0;
894
895 /* presume no rounding */
896 round_inc = 0;
897
898 /* check to see if the first fractional digit is represented */
899 if (exp >= 0 && (size_t)exp < prec)
900 {
901 /* if the digit is 5 or over, round up */
902 if (get_dig(ext_, (size_t)exp) >= 5)
903 round_inc = (is_neg ? -1 : 1);
904 }
905
906 /* stop at the first fractional digit */
907 if (exp <= 0)
908 {
909 /* all digits are fractional */
910 stop_idx = 0;
911 }
912 else if ((size_t)exp < prec)
913 {
914 /* stop at the first fractional digit */
915 stop_idx = (size_t)exp;
916 }
917 else
918 {
919 /* all of the digits are in the whole part */
920 stop_idx = prec;
921 }
922
923 /* convert the integer part digit by digit */
924 if (stop_idx > 0)
925 {
926 /* loop over digits */
927 for (idx = 0 ; ; ++idx)
928 {
929 /* add this digit */
930 if (is_neg)
931 acc -= get_dig(ext_, idx);
932 else
933 acc += get_dig(ext_, idx);
934
935 /* if we have another digit waiting, multiply by ten */
936 if (idx + 1 < stop_idx)
937 {
938 /* make sure this wouldn't overflow */
939 if ((is_neg && acc < long_min_over_10 - round_inc)
940 || (!is_neg && acc > long_max_over_10 - round_inc))
941 {
942 /* the number is too large for an integer */
943 err_throw(VMERR_NUM_OVERFLOW);
944 }
945
946 /* multiply the accumulator by 10 for the next digit */
947 acc *= 10;
948 }
949 else
950 {
951 /* we're done */
952 break;
953 }
954 }
955 }
956
957 /* return the result adjusted for rounding */
958 return acc + round_inc;
959 }
960
961 /* ------------------------------------------------------------------------ */
962 /*
963 * Create a string representation of the number
964 */
cast_to_string(VMG_ vm_obj_id_t self,vm_val_t * new_str) const965 const char *CVmObjBigNum::cast_to_string(VMG_ vm_obj_id_t self,
966 vm_val_t *new_str) const
967 {
968 /* use my static method */
969 return cvt_to_string(vmg_ self, new_str, ext_, 100, -1, -1, -1, 0, 0);
970 }
971
972 /*
973 * convert to a string, storing the result in the given buffer
974 */
cvt_to_string_buf(VMG_ char * buf,size_t buflen,int max_digits,int whole_places,int frac_digits,int exp_digits,ulong flags)975 char *CVmObjBigNum::cvt_to_string_buf(VMG_ char *buf, size_t buflen,
976 int max_digits, int whole_places,
977 int frac_digits, int exp_digits,
978 ulong flags)
979 {
980 /* convert to a string into our buffer */
981 return cvt_to_string_gen(vmg_ 0, ext_, max_digits, whole_places,
982 frac_digits, exp_digits, flags, 0,
983 buf, buflen);
984 }
985
986 /*
987 * Convert to a string, creating a new string object to hold the result
988 */
cvt_to_string(VMG_ vm_obj_id_t self,vm_val_t * new_str,const char * ext,int max_digits,int whole_places,int frac_digits,int exp_digits,ulong flags,vm_val_t * lead_fill)989 const char *CVmObjBigNum::cvt_to_string(VMG_ vm_obj_id_t self,
990 vm_val_t *new_str,
991 const char *ext,
992 int max_digits, int whole_places,
993 int frac_digits, int exp_digits,
994 ulong flags, vm_val_t *lead_fill)
995 {
996 const char *ret;
997
998 /* push a self-reference for gc protection during allocation */
999 G_stk->push()->set_obj(self);
1000
1001 /*
1002 * convert to a string - don't pass in a buffer, since we want to
1003 * create a new string to hold the result
1004 */
1005 ret = cvt_to_string_gen(vmg_ new_str, ext, max_digits, whole_places,
1006 frac_digits, exp_digits, flags, lead_fill, 0, 0);
1007
1008 /* discard our gc protection */
1009 G_stk->discard();
1010
1011 /* return the result */
1012 return ret;
1013 }
1014
1015 /*
1016 * Common routine to create a string representation. If buf is null,
1017 * we'll allocate a new string object, filling in new_str with the
1018 * object reference; otherwise, we'll format into the given buffer.
1019 */
cvt_to_string_gen(VMG_ vm_val_t * new_str,const char * ext,int max_digits,int whole_places,int frac_digits,int exp_digits,ulong flags,vm_val_t * lead_fill,char * buf,size_t buflen)1020 char *CVmObjBigNum::cvt_to_string_gen(VMG_ vm_val_t *new_str,
1021 const char *ext,
1022 int max_digits, int whole_places,
1023 int frac_digits, int exp_digits,
1024 ulong flags, vm_val_t *lead_fill,
1025 char *buf, size_t buflen)
1026 {
1027 int always_sign = ((flags & VMBN_FORMAT_SIGN) != 0);
1028 int always_sign_exp = ((flags & VMBN_FORMAT_EXP_SIGN) != 0);
1029 int always_exp = ((flags & VMBN_FORMAT_EXP) != 0);
1030 int lead_zero = ((flags & VMBN_FORMAT_LEADING_ZERO) != 0);
1031 int always_show_pt = ((flags & VMBN_FORMAT_POINT) != 0);
1032 int exp_caps = ((flags & VMBN_FORMAT_EXP_CAP) != 0);
1033 int pos_lead_space = ((flags & VMBN_FORMAT_POS_SPACE) != 0);
1034 int commas = ((flags & VMBN_FORMAT_COMMAS) != 0);
1035 int euro = ((flags & VMBN_FORMAT_EUROSTYLE) != 0);
1036 size_t req_chars;
1037 int prec = (int)get_prec(ext);
1038 int exp = get_exp(ext);
1039 int dig_before_pt;
1040 int dig_after_pt;
1041 int idx;
1042 unsigned int dig;
1043 char *p;
1044 int exp_carry;
1045 int show_pt;
1046 int whole_padding;
1047 int non_sci_digs;
1048 char decpt_char = (euro ? ',' : '.');
1049 char comma_char = (euro ? '.' : ',');
1050 const char *lead_fill_str = 0;
1051 size_t lead_fill_len;
1052 char *tmp_ext = 0;
1053 uint tmp_hdl;
1054
1055 /* get the fill string, if a value was provided */
1056 if (lead_fill != 0 && lead_fill->typ != VM_NIL)
1057 {
1058 /* get the string value */
1059 lead_fill_str = lead_fill->get_as_string(vmg0_);
1060
1061 /* if it's not a string, it's an error */
1062 if (lead_fill_str == 0)
1063 err_throw(VMERR_STRING_VAL_REQD);
1064
1065 /* read and skip the length prefix */
1066 lead_fill_len = vmb_get_len(lead_fill_str);
1067 lead_fill_str += VMB_LEN;
1068
1069 /* if the length is zero, ignore the lead fill string entirely */
1070 if (lead_fill_len == 0)
1071 lead_fill_str = 0;
1072 }
1073 else
1074 {
1075 /* no lead fill needed */
1076 lead_fill_len = 0;
1077 }
1078
1079 /*
1080 * If we're not required to use exponential notation, but we don't
1081 * have enough max_digits places for the part before the decimal
1082 * point, use exponential anyway. (The number of digits before the
1083 * decimal point is simply the exponent if it's greater than zero,
1084 * or zero otherwise.)
1085 */
1086 if (exp > max_digits)
1087 always_exp = TRUE;
1088
1089 /*
1090 * If we're not required to use exponential notation, but our
1091 * absolute value is so small that we wouldn't show anything
1092 * "0.00000..." (i.e., we'd have too many zeroes after the decimal
1093 * point to show any actual digits of our number), use exponential
1094 * notation. If our exponent is negative, its absolute value is the
1095 * number of zeroes we'd show after the decimal point before the
1096 * first non-zero digit.
1097 */
1098 if (exp < 0
1099 && (-exp > max_digits
1100 || (frac_digits != -1 && -exp > frac_digits)))
1101 always_exp = TRUE;
1102
1103 /* calculate how many digits we'd need in non-scientific notation */
1104 if (exp < 0)
1105 {
1106 /* we have leading zeroes before the first significant digit */
1107 non_sci_digs = -exp + prec;
1108 }
1109 else if (exp > prec)
1110 {
1111 /* we have trailing zeroes after the last significant digit */
1112 non_sci_digs = exp + prec;
1113 }
1114 else
1115 {
1116 /*
1117 * we have no leading or trailing zeroes to represent - only the
1118 * digits actually stored need to be displayed
1119 */
1120 non_sci_digs = prec;
1121 }
1122
1123 /*
1124 * Figure out how much space we need for our string: use the smaller
1125 * of max_digits or the actual space we need for non-scientific
1126 * notation, plus overhead space for the sign, a leading zero, a
1127 * decimal point, an 'E' for the exponent, an exponent sign, and up
1128 * to five digits for the exponent (16-bit integer -> -32768 to
1129 * 32767). Also add one extra digit in case we need to add a digit
1130 * due to rounding.
1131 */
1132 if (max_digits < non_sci_digs)
1133 req_chars = max_digits;
1134 else
1135 req_chars = non_sci_digs;
1136 req_chars += 11;
1137
1138 /*
1139 * Make sure we leave enough room for the lead fill string - if they
1140 * specified a number of whole places, and we're not using
1141 * exponential format, we'll insert lead fill characters before the
1142 * first non-zero whole digit.
1143 */
1144 if (!always_exp && whole_places != -1 && exp < whole_places)
1145 {
1146 int extra;
1147 int char_size;
1148
1149 /*
1150 * if the exponent is negative, we'll pad by the full amount;
1151 * otherwise, we'll just pad by the difference between the
1152 * number of places needed and the exponent
1153 */
1154 extra = whole_places;
1155 if (exp > 0)
1156 extra -= exp;
1157
1158 /*
1159 * Add the extra bytes: one byte per character if we're using
1160 * the default space padding, or up to three bytes per character
1161 * if a lead string was specified (each unicode character can
1162 * take up to three bytes)
1163 */
1164 char_size = (lead_fill_str != 0 ? 3 : 1);
1165 req_chars += extra * char_size;
1166
1167 /*
1168 * add space for each padding character we could insert into a
1169 * comma position (there's at most one comma per three fill
1170 * characters)
1171 */
1172 if (commas)
1173 req_chars += ((extra + 2)/3) * char_size;
1174 }
1175
1176 /*
1177 * If we're using commas, and we're not using scientific notation,
1178 * add space for a comma for each three digits before the decimal
1179 * point
1180 */
1181 if (commas && !always_exp)
1182 {
1183 /* add space for the commas */
1184 req_chars += ((exp + 2) / 3);
1185 }
1186
1187 /*
1188 * if they requested a specific minimum number of exponent digits,
1189 * and that number is greater than the allowance of 5 we already
1190 * provided, add the extra space needed
1191 */
1192 if (exp_digits > 5)
1193 req_chars += (exp_digits - 5);
1194
1195 /*
1196 * If they requested a specific number of digits after the decimal
1197 * point, make sure we have room for those digits.
1198 */
1199 if (frac_digits != -1)
1200 req_chars += frac_digits;
1201
1202 /* check to see if the caller passed in a buffer */
1203 if (buf != 0)
1204 {
1205 /*
1206 * the caller passed in a buffer - if it's not big enough to
1207 * hold the result, return failure
1208 */
1209 if (buflen < req_chars + VMB_LEN)
1210 return 0;
1211 }
1212 else
1213 {
1214 /* no buffer - allocate a new string */
1215 buf = CVmObjString::alloc_str_buf(vmg_ new_str, 0, 0, req_chars);
1216 }
1217
1218 /* check for special values */
1219 switch(get_type(ext))
1220 {
1221 case VMBN_T_NUM:
1222 /* normal number - proceed */
1223 break;
1224
1225 case VMBN_T_NAN:
1226 /* not a number - show "1.#NAN" */
1227 strcpy(buf + VMB_LEN, "1.#NAN");
1228 oswp2(buf, 6);
1229 return buf;
1230
1231 case VMBN_T_INF:
1232 /* positive or negative infinity */
1233 if (get_neg(ext))
1234 {
1235 strcpy(buf + VMB_LEN, "-1.#INF");
1236 oswp2(buf, 7);
1237 }
1238 else
1239 {
1240 strcpy(buf + VMB_LEN, "1.#INF");
1241 oswp2(buf, 6);
1242 }
1243 return buf;
1244
1245 default:
1246 /* other values are not valid */
1247 strcpy(buf + VMB_LEN, "1.#???");
1248 oswp2(buf, 6);
1249 return buf;
1250 }
1251
1252 /*
1253 * Allocate a temporary buffer to contain a copy of the number. At
1254 * most, we'll have to show max_digits of the number, or the current
1255 * precision, whichever is lower.
1256 */
1257 {
1258 int new_prec;
1259
1260 /*
1261 * limit the new precision to the maximum digits to be shown, or
1262 * our existing precision, whichever is lower
1263 */
1264 new_prec = max_digits;
1265 if (prec < new_prec)
1266 new_prec = prec;
1267
1268 /* allocate the space */
1269 alloc_temp_regs(vmg_ (size_t)new_prec, 1, &tmp_ext, &tmp_hdl);
1270
1271 /* copy the value to the temp register, rounding the value */
1272 copy_val(tmp_ext, ext, TRUE);
1273
1274 /* note the new precision */
1275 prec = new_prec;
1276
1277 /* forget the original number and use the rounded version */
1278 ext = tmp_ext;
1279 }
1280
1281 start_over:
1282 /*
1283 * note the exponent, in case we rounded or otherwise adjusted the
1284 * temporary number
1285 */
1286 exp = get_exp(ext);
1287
1288 /* presume we won't carry into the exponent */
1289 exp_carry = FALSE;
1290
1291 /* no whole-part padding yet */
1292 whole_padding = 0;
1293
1294 /*
1295 * Figure out where the decimal point goes in the display. If we're
1296 * displaying in exponential format, we'll always display exactly
1297 * one digit before the decimal point. Otherwise, we'll display the
1298 * number given by our exponent if it's positive, or zero or one
1299 * (depending on lead_zero) if it's negative or zero.
1300 */
1301 if (always_exp)
1302 dig_before_pt = 1;
1303 else if (exp > 0)
1304 dig_before_pt = exp;
1305 else
1306 dig_before_pt = 0;
1307
1308 /*
1309 * if the digits before the decimal point exceed our maximum number
1310 * of digits allowed, we'll need to use exponential format
1311 */
1312 if (dig_before_pt > max_digits)
1313 {
1314 always_exp = TRUE;
1315 goto start_over;
1316 }
1317
1318 /*
1319 * Limit digits after the decimal point according to the maximum
1320 * digits allowed and the number we'll show before the decimal
1321 * point.
1322 */
1323 dig_after_pt = max_digits - dig_before_pt;
1324
1325 /* start writing after the buffer length prefix */
1326 p = buf + VMB_LEN;
1327
1328 /*
1329 * if we're not in exponential mode, add leading spaces for the
1330 * unused whole places
1331 */
1332 if (!always_exp && dig_before_pt < whole_places)
1333 {
1334 int cnt;
1335 size_t src_rem;
1336 utf8_ptr src;
1337 utf8_ptr dst;
1338 int idx;
1339
1340 /* start with the excess whole places */
1341 cnt = whole_places - dig_before_pt;
1342
1343 /* if we're going to add a leading zero, that's one less space */
1344 if (dig_before_pt == 0 && lead_zero)
1345 --cnt;
1346
1347 /*
1348 * Increase the count by the number of comma positions in the
1349 * padding area.
1350 */
1351 if (commas)
1352 {
1353 /* scan over the positions to fill and count commas */
1354 for (idx = dig_before_pt ; idx < whole_places ; ++idx)
1355 {
1356 /* if this is a comma position, note it */
1357 if ((idx % 3) == 0)
1358 ++cnt;
1359 }
1360 }
1361
1362 /* set up our read and write pointers */
1363 src.set((char *)lead_fill_str);
1364 src_rem = lead_fill_len;
1365 dst.set(p);
1366
1367 /* add (and count) the leading spaces */
1368 for ( ; cnt != 0 ; --cnt, ++whole_padding)
1369 {
1370 wchar_t ch;
1371
1372 /*
1373 * if we have a lead fill string, read from it; otherwise,
1374 * just use a space
1375 */
1376 if (lead_fill_str != 0)
1377 {
1378 /* if we've exhausted the string, start over */
1379 if (src_rem == 0)
1380 {
1381 src.set((char *)lead_fill_str);
1382 src_rem = lead_fill_len;
1383 }
1384
1385 /* get the next character */
1386 ch = src.getch();
1387
1388 /* skip this character */
1389 src.inc(&src_rem);
1390 }
1391 else
1392 {
1393 /* no lead fill string - insert a space by default */
1394 ch = ' ';
1395 }
1396
1397 /* write this character */
1398 dst.setch(ch);
1399 }
1400
1401 /* resynchronize from our utf8 pointer */
1402 p = dst.getptr();
1403 }
1404
1405 /*
1406 * if the number is negative, or we're always showing a sign, add
1407 * the sign; if we're not adding a sign, but we're showing a leading
1408 * space for positive numbers, add the leading space
1409 */
1410 if (get_neg(ext))
1411 *p++ = '-';
1412 else if (always_sign)
1413 *p++ = '+';
1414 else if (pos_lead_space)
1415 *p++ = ' ';
1416
1417 /*
1418 * if we have no digits before the decimal, and we're adding a
1419 * leading zero, add it now
1420 */
1421 if (dig_before_pt == 0 && lead_zero)
1422 {
1423 /* add the leading zero */
1424 *p++ = '0';
1425
1426 /* reduce the limit on the digits after the decimal point */
1427 --dig_after_pt;
1428 }
1429
1430 /*
1431 * If we have limited the number of digits that we'll allow after the
1432 * decimal point, due to the limit on the total number of digits and
1433 * the number of digits we need to display before the decimal, start
1434 * over in exponential format to ensure we get the after-decimal
1435 * display we want.
1436 *
1437 * Note that we won't bother going into exponential format if the
1438 * number of digits before the decimal point is zero or one, because
1439 * exponential format won't give us any more room - in such cases we
1440 * simply have an impossible request.
1441 */
1442 if (!always_exp && frac_digits != -1 && dig_after_pt < frac_digits
1443 && dig_before_pt > 1)
1444 {
1445 /* switch to exponential format and start over */
1446 always_exp = TRUE;
1447 goto start_over;
1448 }
1449
1450 /* display the digits before the decimal point */
1451 for (idx = 0 ; idx < dig_before_pt && idx < prec ; ++idx)
1452 {
1453 /*
1454 * if this isn't the first digit, and we're adding commas, and
1455 * an even multiple of three more digits follow, insert a comma
1456 */
1457 if (idx != 0 && commas && !always_exp
1458 && ((dig_before_pt - idx) % 3) == 0)
1459 *p++ = comma_char;
1460
1461 /* get this digit */
1462 dig = get_dig(ext, idx);
1463
1464 /* add it to the string */
1465 *p++ = (dig + '0');
1466 }
1467
1468 /* if we ran out of precision, add zeroes */
1469 for ( ; idx < dig_before_pt ; ++idx)
1470 *p++ = '0';
1471
1472 /*
1473 * Add the decimal point. Show the decimal point unless
1474 * always_show_pt is false, and either frac_digits is zero, or
1475 * frac_digits is -1 and we have no fractional part.
1476 */
1477 if (always_show_pt)
1478 {
1479 /* always showing the point */
1480 show_pt = TRUE;
1481 }
1482 else
1483 {
1484 if (frac_digits > 0)
1485 {
1486 /* we're showing fractional digits - always show a point */
1487 show_pt = TRUE;
1488 }
1489 else if (frac_digits == 0)
1490 {
1491 /* we're showing no fractional digits, so suppress the point */
1492 show_pt = FALSE;
1493 }
1494 else
1495 {
1496 /*
1497 * for now assume we'll show the point; we'll take it back
1498 * out if we don't encounter anything but zeroes
1499 */
1500 show_pt = TRUE;
1501 }
1502 }
1503
1504 /* if we're showing the fractional part, show it */
1505 if (show_pt)
1506 {
1507 int frac_len;
1508 int frac_lim;
1509 char *last_non_zero;
1510
1511 /*
1512 * remember the current position as the last trailing non-zero -
1513 * if we don't find anything but zeroes and decide to remove the
1514 * trailing zeroes, we'll also remove the decimal point by
1515 * coming back here
1516 */
1517 last_non_zero = p;
1518
1519 /* add the point */
1520 *p++ = decpt_char;
1521
1522 /* if we're always showing a decimal point, we can't remove it */
1523 if (always_show_pt)
1524 last_non_zero = p;
1525
1526 /* if frac_digits is -1, there's no limit */
1527 if (frac_digits == -1)
1528 frac_lim = dig_after_pt;
1529 else
1530 frac_lim = frac_digits;
1531
1532 /*
1533 * further limit the fractional digits according to the maximum
1534 * digit allowance
1535 */
1536 if (frac_lim > dig_after_pt)
1537 frac_lim = dig_after_pt;
1538
1539 /* no fractional digits output yet */
1540 frac_len = 0;
1541
1542 /*
1543 * if we haven't yet reached the first non-zero digit, display
1544 * as many zeroes as necessary
1545 */
1546 if (idx == 0 && exp < 0)
1547 {
1548 int cnt;
1549
1550 /*
1551 * display leading zeroes based no the exponent: if exp is
1552 * zero, we don't need any; if exp is -1, we need one; if
1553 * exp is -2, we need two, and so on
1554 */
1555 for (cnt = exp ; cnt != 0 && frac_len < frac_lim ;
1556 ++cnt, ++frac_len)
1557 {
1558 /* add a zero */
1559 *p++ = '0';
1560 }
1561 }
1562
1563 /* add the fractional digits */
1564 for ( ; idx < prec && frac_len < frac_lim ; ++idx, ++frac_len)
1565 {
1566 /* get this digit */
1567 dig = get_dig(ext, idx);
1568
1569 /* add it */
1570 *p++ = (dig + '0');
1571
1572 /*
1573 * if it's not a zero, note the location - if we decide to
1574 * trim trailing zeroes, we'll want to keep at least this
1575 * much, since this is a significant trailing digit
1576 */
1577 if (dig != 0)
1578 last_non_zero = p;
1579 }
1580
1581 /*
1582 * add the trailing zeroes if we ran out of precision before we
1583 * filled the requested number of places
1584 */
1585 if (frac_digits != -1)
1586 {
1587 /* fill out the remaining request length with zeroes */
1588 for ( ; frac_len < frac_lim ; ++frac_len)
1589 *p++ = '0';
1590 }
1591 else
1592 {
1593 char *p1;
1594
1595 /*
1596 * if we're removing the decimal point, we're not showing a
1597 * fractional part after all - so note
1598 */
1599 if (last_non_zero < p && *last_non_zero == decpt_char)
1600 show_pt = FALSE;
1601
1602 /*
1603 * We can use whatever length we like, so remove meaningless
1604 * trailing zeroes. Before we do this, though, make sure we
1605 * aren't rounding up the last trailing zero - if the next
1606 * digit is 5 or higher, we'll round the final zero to a 1.
1607 */
1608 if (p > last_non_zero
1609 && idx < prec
1610 && get_dig(ext, idx) >= 5)
1611 {
1612 /*
1613 * That last zero is significant after all, since we're
1614 * going to round it up to a 1 for display. Don't actually
1615 * do the rounding now; simply note that the last zero is
1616 * significant so that we don't drop the digits leading up
1617 * to it.
1618 */
1619 last_non_zero = p;
1620 }
1621
1622 /*
1623 * We've checked for rounding in the last digit, so we can now
1624 * safely remove meaningless trailing zeroes. If this leaves a
1625 * completely empty buffer, not counting a sign and/or a
1626 * decimal point, it means that we have a fractional number
1627 * that we're showing without an exponent, and the number of
1628 * digits we had for display was insufficient to reach the
1629 * first non-zero fractional digit. In this case, simply show
1630 * '0' (or '0.', if a decimal point is required) as the result.
1631 * To find out, scan for digits.
1632 */
1633 p = last_non_zero;
1634 for (p1 = buf + VMB_LEN ; p1 < p && !is_digit(*p1) ; ++p1) ;
1635
1636 /* if we didn't find any digits, add/insert a '0' */
1637 if (p1 == p)
1638 {
1639 /*
1640 * if there's a decimal point, insert the '0' before it;
1641 * otherwise, just append the zero
1642 */
1643 if (p > buf + VMB_LEN && *(p-1) == decpt_char)
1644 {
1645 *(p-1) = '0';
1646 *p++ = decpt_char;
1647 }
1648 else
1649 *p++ = '0';
1650 }
1651 }
1652 }
1653
1654 /*
1655 * Check for rounding. If another digit remains, and that digit is
1656 * greater than or equal to 5, round up.
1657 */
1658 if (idx < prec && get_dig(ext, idx) >= 5)
1659 {
1660 char *rp;
1661 int need_carry;
1662 int found_pt;
1663 int dig_cnt;
1664
1665 /*
1666 * go back through the number and add one to the last digit,
1667 * carrying as needed
1668 */
1669 for (dig_cnt = 0, found_pt = FALSE, need_carry = TRUE, rp = p - 1 ;
1670 rp >= buf + VMB_LEN ; rp = utf8_ptr::s_dec(rp))
1671 {
1672 /* if this is a digit, bump it up */
1673 if (is_digit(*rp))
1674 {
1675 /* count it */
1676 ++dig_cnt;
1677
1678 /* if it's 9, we'll have to carry; otherwise it's easy */
1679 if (*rp == '9')
1680 {
1681 /* set it to zero and keep going to do the carry */
1682 *rp = '0';
1683
1684 /*
1685 * if we haven't found the decimal point yet, and
1686 * we're not required to show a certain number of
1687 * fractional digits, we can simply remove this
1688 * trailing zero
1689 */
1690 if (show_pt && !found_pt && frac_digits == -1)
1691 {
1692 /* remove the trailing zero */
1693 p = utf8_ptr::s_dec(p);
1694
1695 /* remove it from the digit count */
1696 --dig_cnt;
1697 }
1698 }
1699 else
1700 {
1701 /* bump it up one */
1702 ++(*rp);
1703
1704 /* no more carrying is needed */
1705 need_carry = FALSE;
1706
1707 /* we don't need to look any further */
1708 break;
1709 }
1710 }
1711 else if (*rp == decpt_char)
1712 {
1713 /* note that we found the decimal point */
1714 found_pt = TRUE;
1715 }
1716 }
1717
1718 /*
1719 * If we got to the start of the number and we still need a
1720 * carry, we must have had all 9's. In this case, we need to
1721 * redo the entire number, since all of the layout (commas,
1722 * leading spaces, etc) can change - it's way too much work to
1723 * try to back-patch all of this stuff, so we'll just bump the
1724 * number up and reformat it from scratch.
1725 */
1726 if (need_carry)
1727 {
1728 int carry;
1729
1730 /*
1731 * clear the digit that provoked the rounding - we don't
1732 * want to round again on the next iteration
1733 */
1734 set_dig(tmp_ext, idx, 0);
1735
1736 /* round up the number starting at the previous digit */
1737 for (carry = TRUE ; idx != 0 ; )
1738 {
1739 /* move to the previous digit */
1740 --idx;
1741
1742 /* if this digit is a 9, we'll need to carry */
1743 if (get_dig(tmp_ext, idx) == 9)
1744 {
1745 /* adjust this digit and keep going */
1746 set_dig(tmp_ext, idx, 0);
1747 }
1748 else
1749 {
1750 /* bump this digit up one */
1751 set_dig(tmp_ext, idx, get_dig(tmp_ext, idx) + 1);
1752
1753 /* we're done */
1754 carry = FALSE;
1755 break;
1756 }
1757 }
1758
1759 /* if we need to carry one more place, shift it */
1760 if (carry)
1761 {
1762 /* shift the number */
1763 shift_right(tmp_ext, 1);
1764
1765 /* adjust the exponent accordingly */
1766 set_exp(tmp_ext, get_exp(tmp_ext) + 1);
1767
1768 /* insert the leading 1 */
1769 set_dig(tmp_ext, 0, 1);
1770 }
1771
1772 /*
1773 * if this pushes us over the maximum digit range, switch to
1774 * scientific notation
1775 */
1776 if (dig_cnt + 1 > max_digits)
1777 always_exp = TRUE;
1778
1779 /* go back and start over */
1780 goto start_over;
1781 }
1782 }
1783
1784 /* add the exponent */
1785 if (always_exp)
1786 {
1787 int disp_exp;
1788
1789 /* add the 'E' */
1790 *p++ = (exp_caps ? 'E' : 'e');
1791
1792 /*
1793 * calculate the display exponent - it's one less than the
1794 * actual exponent, because we put the point after one digit
1795 */
1796 disp_exp = exp - 1;
1797
1798 /*
1799 * if we carried into the exponent due to rounding, increase the
1800 * exponent by one
1801 */
1802 if (exp_carry)
1803 ++disp_exp;
1804
1805 /* add the sign */
1806 if (disp_exp < 0)
1807 {
1808 *p++ = '-';
1809 disp_exp = -disp_exp;
1810 }
1811 else if (always_sign_exp)
1812 *p++ = '+';
1813
1814 /*
1815 * if the exponent is zero, just put zero (unless a more
1816 * specific format was requested)
1817 */
1818 if (disp_exp == 0 && exp_digits == -1)
1819 {
1820 /* add the zero exponent */
1821 *p++ = '0';
1822 }
1823 else
1824 {
1825 char buf[20];
1826 char *ep;
1827 int dig_cnt;
1828
1829 /* build the exponent in reverse */
1830 for (dig_cnt = 0, ep = buf + sizeof(buf) ; disp_exp != 0 ;
1831 disp_exp /= 10, ++dig_cnt)
1832 {
1833 /* move back one character */
1834 --ep;
1835
1836 /* add one digit */
1837 *ep = (disp_exp % 10) + '0';
1838 }
1839
1840 /* if necessary, add leading zeroes to the exponent */
1841 if (exp_digits != -1 && exp_digits > dig_cnt)
1842 {
1843 for ( ; dig_cnt < exp_digits ; ++dig_cnt)
1844 *p++ = '0';
1845 }
1846
1847 /* copy the exponent into the output */
1848 for ( ; ep < buf + sizeof(buf) ; ++ep)
1849 *p++ = *ep;
1850 }
1851 }
1852
1853 /* set the string length */
1854 vmb_put_len(buf, p - (buf + VMB_LEN));
1855
1856 /* if we allocated a temporary register, free it */
1857 if (tmp_ext != 0)
1858 release_temp_regs(vmg_ 1, tmp_hdl);
1859
1860 /* return the string buffer */
1861 return buf;
1862 }
1863
1864 /* ------------------------------------------------------------------------ */
1865 /*
1866 * Shift the value left (multiply by 10^shift)
1867 */
shift_left(char * ext,unsigned int shift)1868 void CVmObjBigNum::shift_left(char *ext, unsigned int shift)
1869 {
1870 size_t prec = get_prec(ext);
1871 size_t i;
1872
1873 /* do nothing with a zero shift */
1874 if (shift == 0)
1875 return;
1876
1877 /* if it's an even shift, it's especially easy */
1878 if ((shift & 1) == 0)
1879 {
1880 /* simply move the bytes left by the required amount */
1881 for (i = 0 ; i + shift/2 < (prec+1)/2 ; ++i)
1882 ext[VMBN_MANT + i] = ext[VMBN_MANT + i + shift/2];
1883
1884 /* zero the remaining digits */
1885 for ( ; i < (prec+1)/2 ; ++i)
1886 ext[VMBN_MANT + i] = 0;
1887
1888 /*
1889 * be sure to zero the last digit - if we have an odd precision,
1890 * we will have copied the garbage digit from the final
1891 * half-byte
1892 */
1893 set_dig(ext, prec - shift, 0);
1894 }
1895 else
1896 {
1897 /* apply the shift to each digit */
1898 for (i = 0 ; i + shift < prec ; ++i)
1899 {
1900 unsigned int dig;
1901
1902 /* get this source digit */
1903 dig = get_dig(ext, i + shift);
1904
1905 /* set this destination digit */
1906 set_dig(ext, i, dig);
1907 }
1908
1909 /* zero the remaining digits */
1910 for ( ; i < prec ; ++i)
1911 set_dig(ext, i, 0);
1912 }
1913 }
1914
1915 /*
1916 * Shift the value right (divide by 10^shift)
1917 */
shift_right(char * ext,unsigned int shift)1918 void CVmObjBigNum::shift_right(char *ext, unsigned int shift)
1919 {
1920 size_t prec = get_prec(ext);
1921 size_t i;
1922
1923 /* if it's an even shift, it's especially easy */
1924 if ((shift & 1) == 0)
1925 {
1926 /* simply move the bytes left by the required amount */
1927 for (i = (prec+1)/2 ; i > shift/2 ; --i)
1928 ext[VMBN_MANT + i-1] = ext[VMBN_MANT + i-1 - shift/2];
1929
1930 /* zero the leading digits */
1931 for ( ; i > 0 ; --i)
1932 ext[VMBN_MANT + i-1] = 0;
1933 }
1934 else
1935 {
1936 /* apply the shift to each digit */
1937 for (i = prec ; i > shift ; --i)
1938 {
1939 unsigned int dig;
1940
1941 /* get this source digit */
1942 dig = get_dig(ext, i-1 - shift);
1943
1944 /* set this destination digit */
1945 set_dig(ext, i-1, dig);
1946 }
1947
1948 /* zero the remaining digits */
1949 for ( ; i >0 ; --i)
1950 set_dig(ext, i-1, 0);
1951 }
1952 }
1953
1954 /*
1955 * Increment a number
1956 */
increment_abs(char * ext)1957 void CVmObjBigNum::increment_abs(char *ext)
1958 {
1959 size_t idx;
1960 int exp = get_exp(ext);
1961 int carry;
1962
1963 /* start at the one's place, if represented */
1964 idx = (exp <= 0 ? 0 : (size_t)exp);
1965
1966 /*
1967 * if the units digit is to the right of the number (i.e., the
1968 * number's scale is large), there's nothing to do
1969 */
1970 if (idx > get_prec(ext))
1971 return;
1972
1973 /* increment digits */
1974 for (carry = TRUE ; idx != 0 ; )
1975 {
1976 int dig;
1977
1978 /* move to the next digit */
1979 --idx;
1980
1981 /* get the digit value */
1982 dig = get_dig(ext, idx);
1983
1984 /* increment it, checking for carry */
1985 if (dig == 9)
1986 {
1987 /* increment it to zero and keep going to carry */
1988 set_dig(ext, idx, 0);
1989 }
1990 else
1991 {
1992 /* increment this digit */
1993 set_dig(ext, idx, dig + 1);
1994
1995 /* there's no carry out */
1996 carry = FALSE;
1997
1998 /* done */
1999 break;
2000 }
2001 }
2002
2003 /* if we carried past the end of the number, insert the leading 1 */
2004 if (carry)
2005 {
2006 /*
2007 * if we still haven't reached the units position, shift right
2008 * until we do
2009 */
2010 while (get_exp(ext) < 0)
2011 {
2012 /* shift it right and adjust the exponent */
2013 shift_right(ext, 1);
2014 set_exp(ext, get_exp(ext) + 1);
2015 }
2016
2017 /* shift the number right and adjust the exponent */
2018 shift_right(ext, 1);
2019 set_exp(ext, get_exp(ext) + 1);
2020
2021 /* insert the leading 1 */
2022 set_dig(ext, 0, 1);
2023 }
2024
2025 /* we know the value is now non-zero */
2026 ext[VMBN_FLAGS] &= ~VMBN_F_ZERO;
2027 }
2028
2029 /*
2030 * Determine if the fractional part is zero
2031 */
is_frac_zero(const char * ext)2032 int CVmObjBigNum::is_frac_zero(const char *ext)
2033 {
2034 size_t idx;
2035 int exp = get_exp(ext);
2036 size_t prec = get_prec(ext);
2037
2038 /* start at the first fractional digit, if represented */
2039 idx = (exp <= 0 ? 0 : (size_t)exp);
2040
2041 /* scan the digits for a non-zero digit */
2042 for ( ; idx < prec ; ++idx)
2043 {
2044 /* if this digit is non-zero, the fraction is non-zero */
2045 if (get_dig(ext, idx) != 0)
2046 return FALSE;
2047 }
2048
2049 /* we didn't find any non-zero fractional digits */
2050 return TRUE;
2051 }
2052
2053 /*
2054 * Normalize a number - shift it so that the first digit is non-zero.
2055 * If the number is zero, set the exponent to zero and clear the sign
2056 * bit.
2057 */
normalize(char * ext)2058 void CVmObjBigNum::normalize(char *ext)
2059 {
2060 int idx;
2061 int prec = get_prec(ext);
2062 int all_zero;
2063 int nonzero_idx;
2064
2065 /* no work is necessary for anything but ordinary numbers */
2066 if (is_nan(ext))
2067 return;
2068
2069 /* check for an all-zero mantissa */
2070 for (all_zero = TRUE, idx = 0 ; idx < prec ; ++idx)
2071 {
2072 /* check this digit */
2073 if (get_dig(ext, idx) != 0)
2074 {
2075 /* note the location of the first non-zero digit */
2076 nonzero_idx = idx;
2077
2078 /* note that the number isn't all zeroes */
2079 all_zero = FALSE;
2080
2081 /* no need to keep searching */
2082 break;
2083 }
2084 }
2085
2086 /* if it's zero, set the canonical zero format */
2087 if (all_zero)
2088 {
2089 /* set the value to zero */
2090 set_zero(ext);
2091 }
2092 else
2093 {
2094 /* clear the zero flag */
2095 ext[VMBN_FLAGS] &= ~VMBN_F_ZERO;
2096
2097 /* shift the mantissa left to make the first digit non-zero */
2098 if (nonzero_idx != 0)
2099 shift_left(ext, nonzero_idx);
2100
2101 /* decrease the exponent to account for the mantissa shift */
2102 set_exp(ext, get_exp(ext) - nonzero_idx);
2103 }
2104 }
2105
2106 /*
2107 * Round the value up - increments the least significant digit
2108 */
round_up_abs(char * ext)2109 void CVmObjBigNum::round_up_abs(char *ext)
2110 {
2111 int idx;
2112 int carry;
2113
2114 /*
2115 * Scan from least significant and apply the rounding. Keep going
2116 * until we reach the most significant digit.
2117 */
2118 for (carry = TRUE, idx = get_prec(ext) ; idx != 0 ; )
2119 {
2120 int dig;
2121
2122 /* move to the next position */
2123 --idx;
2124
2125 /* get the digit at this position */
2126 dig = get_dig(ext, idx);
2127
2128 /* check to see if we'll need to carry past this digit */
2129 if (dig == 9)
2130 {
2131 /* set it to zero and keep going to do the carry */
2132 set_dig(ext, idx, 0);
2133 }
2134 else
2135 {
2136 /* increment this digit */
2137 set_dig(ext, idx, dig + 1);
2138
2139 /* no need to carry - note it and stop looping */
2140 carry = FALSE;
2141 break;
2142 }
2143 }
2144
2145 /*
2146 * if we carried past the most significant digit, we must shift the
2147 * value right, dropping the least significant digit, and insert a
2148 * leading 1
2149 */
2150 if (carry)
2151 {
2152 /* shift the mantissa */
2153 shift_right(ext, 1);
2154
2155 /* compensate for the shift in the exponent */
2156 set_exp(ext, get_exp(ext) + 1);
2157
2158 /* insert the leading 1 */
2159 set_dig(ext, 0, 1);
2160 }
2161
2162 /* we know the value is non-zero now */
2163 ext[VMBN_FLAGS] &= ~VMBN_F_ZERO;
2164 }
2165
2166 /*
2167 * Copy a value, extending with zeroes if expanding, or truncating or
2168 * rounding, as desired, if the precision changes
2169 */
copy_val(char * dst,const char * src,int round)2170 void CVmObjBigNum::copy_val(char *dst, const char *src, int round)
2171 {
2172 size_t src_prec = get_prec(src);
2173 size_t dst_prec = get_prec(dst);
2174
2175 /* check to see if we're growing or shrinking */
2176 if (dst_prec > src_prec)
2177 {
2178 /*
2179 * it's growing - copy the entire old mantissa, plus the flags,
2180 * sign, and exponent
2181 */
2182 memcpy(dst + VMBN_EXP, src + VMBN_EXP,
2183 (VMBN_MANT - VMBN_EXP) + (src_prec + 1)/2);
2184
2185 /* clear the balance of the mantissa */
2186 memset(dst + VMBN_MANT + (src_prec + 1)/2,
2187 0, (dst_prec + 1)/2 - (src_prec + 1)/2);
2188
2189 /*
2190 * clear the digit just after the last digit we copied - we
2191 * might have missed this with the memset, since it only deals
2192 * with even-numbered pairs of digits
2193 */
2194 set_dig(dst, src_prec, 0);
2195 }
2196 else
2197 {
2198 /* it's shrinking - truncate the mantissa to the new length */
2199 memcpy(dst + VMBN_EXP, src + VMBN_EXP,
2200 (VMBN_MANT - VMBN_EXP) + (dst_prec + 1)/2);
2201
2202 /* check for rounding */
2203 if (round && dst_prec < src_prec && get_dig(src, dst_prec) >= 5)
2204 {
2205 /* round the value */
2206 round_up_abs(dst);
2207 }
2208 }
2209 }
2210
2211 /*
2212 * Multiply by an integer constant value
2213 */
mul_by_long(char * ext,unsigned long val)2214 void CVmObjBigNum::mul_by_long(char *ext, unsigned long val)
2215 {
2216 size_t idx;
2217 size_t prec = get_prec(ext);
2218 unsigned long carry = 0;
2219 int trail_dig = 0;
2220
2221 /*
2222 * start at the least significant digit and work up through the
2223 * digits
2224 */
2225 for (idx = prec ; idx != 0 ; )
2226 {
2227 unsigned long prod;
2228
2229 /* move to the next digit */
2230 --idx;
2231
2232 /*
2233 * compute the product of this digit and the given value, and
2234 * add in the carry from the last digit
2235 */
2236 prod = (val * get_dig(ext, idx)) + carry;
2237
2238 /* set this digit to the residue mod 10 */
2239 set_dig(ext, idx, prod % 10);
2240
2241 /* carry the rest */
2242 carry = prod / 10;
2243 }
2244
2245 /* if we have a carry left over, shift it in */
2246 while (carry != 0)
2247 {
2248 /* remember the digit we're dropping */
2249 trail_dig = get_dig(ext, prec - 1);
2250
2251 /* shift the number and adjust the exponent */
2252 shift_right(ext, 1);
2253 set_exp(ext, get_exp(ext) + 1);
2254
2255 /* shift in this digit of the carry */
2256 set_dig(ext, 0, carry % 10);
2257
2258 /* take it out of the carry */
2259 carry /= 10;
2260 }
2261
2262 /* round up if the dropped trailing digit is 5 or higher */
2263 if (trail_dig >= 5)
2264 round_up_abs(ext);
2265
2266 /* normalize the result */
2267 normalize(ext);
2268 }
2269
2270 /*
2271 * Divide by an integer constant value
2272 */
div_by_long(char * ext,unsigned long val)2273 void CVmObjBigNum::div_by_long(char *ext, unsigned long val)
2274 {
2275 size_t in_idx;
2276 size_t out_idx;
2277 int sig;
2278 size_t prec = get_prec(ext);
2279 unsigned long rem = 0;
2280
2281 /*
2282 * start at the most significant digit and work down
2283 */
2284 for (rem = 0, sig = FALSE, in_idx = out_idx = 0 ;
2285 in_idx < prec || out_idx < prec ; ++in_idx)
2286 {
2287 long quo;
2288
2289 /*
2290 * shift this digit into the remainder - if we're past the end
2291 * of the number, shift in an implied trailing zero
2292 */
2293 rem *= 10;
2294 rem += (in_idx < prec ? get_dig(ext, in_idx) : 0);
2295
2296 /* calculate the quotient */
2297 quo = rem / val;
2298
2299 /* if the quotient is non-zero, we've found a significant digit */
2300 if (quo != 0)
2301 sig = TRUE;
2302
2303 /*
2304 * if we've found a significant digit, store it; otherwise, just
2305 * reduce the exponent to account for an implied leading zero
2306 * that we won't actually store
2307 */
2308 if (sig)
2309 {
2310 /* store the digit */
2311 set_dig(ext, out_idx, (int)quo);
2312
2313 /* move on to the next output digit */
2314 ++out_idx;
2315 }
2316 else
2317 {
2318 /* all leading zeroes so far - adjust the exponent */
2319 set_exp(ext, get_exp(ext) - 1);
2320 }
2321
2322 /* calculate the remainder */
2323 rem %= val;
2324
2325 /*
2326 * if we've reached the last input digit and the remainder is
2327 * zero, we're done - fill out the rest of the number with
2328 * trailing zeroes and stop looping
2329 */
2330 if (rem == 0 && in_idx >= prec)
2331 {
2332 /* check to see if we have any significant digits */
2333 if (sig)
2334 {
2335 /* fill out the rest of the number with zeroes */
2336 for ( ; out_idx < prec ; ++out_idx)
2337 set_dig(ext, out_idx, 0);
2338 }
2339 else
2340 {
2341 /* no significant digits - the result is zero */
2342 set_zero(ext);
2343 }
2344
2345 /* we have our result */
2346 break;
2347 }
2348 }
2349
2350 /*
2351 * Round up if the next digit that we can't store is 5 or higher.
2352 * The next digit can be calculated by shifting in the implied
2353 * trailing zero (i.e., multiplying the remainder by 10 and adding
2354 * zero) then dividing it by the divisor.
2355 */
2356 if ((rem * 10)/val >= 5)
2357 round_up_abs(ext);
2358
2359 /* normalize the result */
2360 normalize(ext);
2361 }
2362
2363 /* ------------------------------------------------------------------------ */
2364 /*
2365 * save to a file
2366 */
save_to_file(VMG_ class CVmFile * fp)2367 void CVmObjBigNum::save_to_file(VMG_ class CVmFile *fp)
2368 {
2369 size_t prec;
2370
2371 /* get our precision */
2372 prec = get_prec(ext_);
2373
2374 /* write the data */
2375 fp->write_bytes(ext_, calc_alloc(prec));
2376 }
2377
2378 /*
2379 * restore from a file
2380 */
restore_from_file(VMG_ vm_obj_id_t,CVmFile * fp,CVmObjFixup *)2381 void CVmObjBigNum::restore_from_file(VMG_ vm_obj_id_t,
2382 CVmFile *fp, CVmObjFixup *)
2383 {
2384 size_t prec;
2385
2386 /* read the precision */
2387 prec = fp->read_uint2();
2388
2389 /* free any existing extension */
2390 if (ext_ != 0)
2391 {
2392 G_mem->get_var_heap()->free_mem(ext_);
2393 ext_ = 0;
2394 }
2395
2396 /* allocate the space */
2397 alloc_bignum(vmg_ prec);
2398
2399 /* store our precision */
2400 set_prec(ext_, prec);
2401
2402 /* read the contents */
2403 fp->read_bytes(ext_ + VMBN_EXP, calc_alloc(prec) - VMBN_EXP);
2404 }
2405
2406
2407 /* ------------------------------------------------------------------------ */
2408 /*
2409 * set a property
2410 */
set_prop(VMG_ class CVmUndo *,vm_obj_id_t,vm_prop_id_t,const vm_val_t *)2411 void CVmObjBigNum::set_prop(VMG_ class CVmUndo *,
2412 vm_obj_id_t, vm_prop_id_t,
2413 const vm_val_t *)
2414 {
2415 /* we have no properties to set */
2416 err_throw(VMERR_INVALID_SETPROP);
2417 }
2418
2419 /*
2420 * get a static property
2421 */
call_stat_prop(VMG_ vm_val_t * result,const uchar ** pc_ptr,uint * argc,vm_prop_id_t prop)2422 int CVmObjBigNum::call_stat_prop(VMG_ vm_val_t *result, const uchar **pc_ptr,
2423 uint *argc, vm_prop_id_t prop)
2424 {
2425 /* translate the property into a function vector index */
2426 switch(G_meta_table
2427 ->prop_to_vector_idx(metaclass_reg_->get_reg_idx(), prop))
2428 {
2429 case VMOBJBN_GET_PI:
2430 return s_getp_pi(vmg_ result, argc);
2431
2432 case VMOBJBN_GET_E:
2433 return s_getp_e(vmg_ result, argc);
2434
2435 default:
2436 /*
2437 * we don't define this one - inherit the default from the base
2438 * object metaclass
2439 */
2440 return CVmObject::call_stat_prop(vmg_ result, pc_ptr, argc, prop);
2441 }
2442 }
2443
2444 /*
2445 * get a property
2446 */
get_prop(VMG_ vm_prop_id_t prop,vm_val_t * val,vm_obj_id_t self,vm_obj_id_t * source_obj,uint * argc)2447 int CVmObjBigNum::get_prop(VMG_ vm_prop_id_t prop, vm_val_t *val,
2448 vm_obj_id_t self, vm_obj_id_t *source_obj,
2449 uint *argc)
2450 {
2451 unsigned short func_idx;
2452
2453 /* translate the property into a function vector index */
2454 func_idx = G_meta_table
2455 ->prop_to_vector_idx(metaclass_reg_->get_reg_idx(), prop);
2456
2457 /* call the function */
2458 if ((this->*func_table_[func_idx])(vmg_ self, val, argc))
2459 {
2460 *source_obj = metaclass_reg_->get_class_obj(vmg0_);
2461 return TRUE;
2462 }
2463
2464 /* inherit default handling */
2465 return CVmObject::get_prop(vmg_ prop, val, self, source_obj, argc);
2466 }
2467
2468 /*
2469 * Property evaluator - formatString
2470 */
getp_format(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)2471 int CVmObjBigNum::getp_format(VMG_ vm_obj_id_t self,
2472 vm_val_t *retval, uint *argc)
2473 {
2474 int orig_argc = (argc != 0 ? *argc : 0);
2475 int max_digits;
2476 uint flags = 0;
2477 int whole_places = -1;
2478 int frac_digits = -1;
2479 int exp_digits = -1;
2480 vm_val_t *lead_fill = 0;
2481 static CVmNativeCodeDesc desc(1, 5);
2482
2483 /* check arguments */
2484 if (get_prop_check_argc(retval, argc, &desc))
2485 return TRUE;
2486
2487 /* get the maximum digit count */
2488 max_digits = CVmBif::pop_int_val(vmg0_);
2489
2490 /* get the flags */
2491 if (orig_argc >= 2)
2492 flags = CVmBif::pop_int_val(vmg0_);
2493
2494 /* get the whole places */
2495 if (orig_argc >= 3)
2496 whole_places = CVmBif::pop_int_val(vmg0_);
2497
2498 /* get the fraction digits */
2499 if (orig_argc >= 4)
2500 frac_digits = CVmBif::pop_int_val(vmg0_);
2501
2502 /* get the exponent digits */
2503 if (orig_argc >= 5)
2504 exp_digits = CVmBif::pop_int_val(vmg0_);
2505
2506 /*
2507 * get the lead fill string if provided (leave it on the stack to
2508 * protect against garbage collection)
2509 */
2510 if (orig_argc >= 6)
2511 lead_fill = G_stk->get(0);
2512
2513 /* format the number, which will build the return string */
2514 cvt_to_string(vmg_ self, retval, ext_, max_digits, whole_places,
2515 frac_digits, exp_digits, flags, lead_fill);
2516
2517 /* discard the lead fill string, if provided */
2518 if (lead_fill != 0)
2519 G_stk->discard();
2520
2521 /* handled */
2522 return TRUE;
2523 }
2524
2525 /*
2526 * Property eval - equal after rounding
2527 */
getp_equal_rnd(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)2528 int CVmObjBigNum::getp_equal_rnd(VMG_ vm_obj_id_t self,
2529 vm_val_t *retval, uint *argc)
2530 {
2531 vm_val_t val2;
2532 static CVmNativeCodeDesc desc(1);
2533
2534 /* check arguments */
2535 if (get_prop_check_argc(retval, argc, &desc))
2536 return TRUE;
2537
2538 /* pop the value to compare */
2539 G_stk->pop(&val2);
2540
2541 /* convert it to BigNumber */
2542 if (!cvt_to_bignum(vmg_ self, &val2))
2543 {
2544 /* it's not a BigNumber, so it's not equal */
2545 retval->set_nil();
2546 }
2547 else
2548 {
2549 int ret;
2550
2551 /* test for equality */
2552 ret = compute_eq_round(vmg_ ext_, get_objid_ext(vmg_ val2.val.obj));
2553
2554 /* set the return value */
2555 retval->set_logical(ret);
2556 }
2557
2558 /* handled */
2559 return TRUE;
2560 }
2561
2562 /*
2563 * property eval - get the precision
2564 */
getp_get_prec(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)2565 int CVmObjBigNum::getp_get_prec(VMG_ vm_obj_id_t self,
2566 vm_val_t *retval, uint *argc)
2567 {
2568 static CVmNativeCodeDesc desc(0);
2569
2570 /* check arguments */
2571 if (get_prop_check_argc(retval, argc, &desc))
2572 return TRUE;
2573
2574 /* return an integer giving my precision */
2575 retval->set_int(get_prec(ext_));
2576
2577 /* handled */
2578 return TRUE;
2579 }
2580
2581 /*
2582 * property eval - set the precision
2583 */
getp_set_prec(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)2584 int CVmObjBigNum::getp_set_prec(VMG_ vm_obj_id_t self,
2585 vm_val_t *retval, uint *argc)
2586 {
2587 size_t digits;
2588 static CVmNativeCodeDesc desc(1);
2589
2590 /* check arguments */
2591 if (get_prop_check_argc(retval, argc, &desc))
2592 return TRUE;
2593
2594 /* get the number of digits for rounding */
2595 digits = CVmBif::pop_int_val(vmg0_);
2596
2597 /* if I'm not a number, return myself unchanged */
2598 if (is_nan(ext_))
2599 {
2600 retval->set_obj(self);
2601 return TRUE;
2602 }
2603
2604 /* push a self-reference while we're working */
2605 G_stk->push()->set_obj(self);
2606
2607 /* create the rounded value */
2608 round_val(vmg_ retval, ext_, digits, TRUE);
2609
2610 /* remove my self-reference */
2611 G_stk->discard();
2612
2613 /* handled */
2614 return TRUE;
2615 }
2616
2617 /*
2618 * get pi (static method)
2619 */
s_getp_pi(VMG_ vm_val_t * val,uint * argc)2620 int CVmObjBigNum::s_getp_pi(VMG_ vm_val_t *val, uint *argc)
2621 {
2622 size_t prec;
2623 char *new_ext;
2624 const char *pi;
2625 static CVmNativeCodeDesc desc(1);
2626
2627 /* check arguments */
2628 if (get_prop_check_argc(val, argc, &desc))
2629 return TRUE;
2630
2631 /* get the precision argument */
2632 prec = CVmBif::pop_int_val(vmg0_);
2633
2634 /* allocate the result */
2635 val->set_obj(create(vmg_ FALSE, prec));
2636 new_ext = get_objid_ext(vmg_ val->val.obj);
2637
2638 /* cache pi to the required precision */
2639 pi = cache_pi(vmg_ prec);
2640
2641 /* return the value */
2642 copy_val(new_ext, pi, TRUE);
2643
2644 /* handled */
2645 return TRUE;
2646 }
2647
2648 /*
2649 * get e (static method)
2650 */
s_getp_e(VMG_ vm_val_t * val,uint * argc)2651 int CVmObjBigNum::s_getp_e(VMG_ vm_val_t *val, uint *argc)
2652 {
2653 size_t prec;
2654 char *new_ext;
2655 const char *e;
2656 static CVmNativeCodeDesc desc(1);
2657
2658 /* check arguments */
2659 if (get_prop_check_argc(val, argc, &desc))
2660 return TRUE;
2661
2662 /* get the precision argument */
2663 prec = CVmBif::pop_int_val(vmg0_);
2664
2665 /* allocate the result */
2666 val->set_obj(create(vmg_ FALSE, prec));
2667 new_ext = get_objid_ext(vmg_ val->val.obj);
2668
2669 /* cache e to the required precision */
2670 e = cache_e(vmg_ prec);
2671
2672 /* return the value */
2673 copy_val(new_ext, e, TRUE);
2674
2675 /* handled */
2676 return TRUE;
2677 }
2678
2679 /*
2680 * Set up for a zero-argument operation that returns a BigNumber result.
2681 * Returns true if the argument check indicates that the caller should
2682 * simply return to its caller, false if the caller should proceed.
2683 *
2684 * After checking the argument count, we'll proceed to set up the return
2685 * value as per setup_getp_retval().
2686 */
setup_getp_0(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc,char ** new_ext)2687 int CVmObjBigNum::setup_getp_0(VMG_ vm_obj_id_t self, vm_val_t *retval,
2688 uint *argc, char **new_ext)
2689 {
2690 static CVmNativeCodeDesc desc(0);
2691
2692 /* check arguments */
2693 if (get_prop_check_argc(retval, argc, &desc))
2694 return TRUE;
2695
2696 /* set up the return value */
2697 return setup_getp_retval(vmg_ self, retval, new_ext, get_prec(ext_));
2698 }
2699
2700 /*
2701 * Set up for a one-argument operation that takes a BigNumber value as
2702 * the argument and returns a BigNumber result. Does the work of
2703 * setup_getp_0, but also pops the argument value and converts it to a
2704 * BigNumber (throwing an error if the value is not convertible).
2705 *
2706 * Fills in val2 with the argument value, and fills in *ext2 with val2's
2707 * extension buffer.
2708 *
2709 * The result value will have the larger of the precisions of self and
2710 * the other value, unless use_self_prec is set, in which case we'll use
2711 * self's precision unconditionally.
2712 *
2713 * If either argument is not a number, we'll set the return value to the
2714 * not-a-number argument unchanged, and return true.
2715 */
setup_getp_1(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc,char ** new_ext,vm_val_t * val2,const char ** ext2,int use_self_prec)2716 int CVmObjBigNum::setup_getp_1(VMG_ vm_obj_id_t self,
2717 vm_val_t *retval, uint *argc,
2718 char **new_ext,
2719 vm_val_t *val2, const char **ext2,
2720 int use_self_prec)
2721 {
2722 size_t prec = get_prec(ext_);
2723 static CVmNativeCodeDesc desc(1);
2724
2725 /* check arguments */
2726 if (get_prop_check_argc(retval, argc, &desc))
2727 return TRUE;
2728
2729 /* pop the argument value */
2730 G_stk->pop(val2);
2731
2732 /* convert it to BigNumber */
2733 if (!cvt_to_bignum(vmg_ self, val2))
2734 {
2735 /* it's not a BigNumber - throw an error */
2736 err_throw(VMERR_BAD_TYPE_BIF);
2737 }
2738
2739 /* get the other value's extension */
2740 *ext2 = get_objid_ext(vmg_ val2->val.obj);
2741
2742 /* if the other value is not a number, return it as the result */
2743 if (is_nan(*ext2))
2744 {
2745 retval->set_obj(val2->val.obj);
2746 return TRUE;
2747 }
2748
2749 /*
2750 * if val2's precision is higher than ours, use it, unless we've
2751 * been specifically told to use our own precision for the result
2752 */
2753 if (!use_self_prec && get_prec(*ext2) > prec)
2754 prec = get_prec(*ext2);
2755
2756 /* push the other result to protect it from garbage collection */
2757 G_stk->push(val2);
2758
2759 /* allocate the return value */
2760 if (setup_getp_retval(vmg_ self, retval, new_ext, prec))
2761 {
2762 /* discard the val2 reference we pushed for gc protection */
2763 G_stk->discard();
2764
2765 /* tell the caller we're done */
2766 return TRUE;
2767 }
2768
2769 /* tell the caller to proceed */
2770 return FALSE;
2771 }
2772
2773 /*
2774 * Set up for an operation that returns a BigNumber result. Returns
2775 * true if we finished the operation, in which case the caller should
2776 * simply return; returns false if the operation should still be carried
2777 * out, in which case the caller should proceed as normal.
2778 *
2779 * If the 'self' value is not-a-number, we'll return it as the result
2780 * (and return true to indicate that no further processing is required).
2781 *
2782 * If we return false, we'll have pushed a reference to 'self' onto the
2783 * stack for protection against garbage collection. The caller must
2784 * discard this reference before returning. We push no such
2785 * self-reference if we return true.
2786 *
2787 * In addition, if we return false, we'll fill in '*new_ext' with a
2788 * pointer to the extension buffer for the return value object that we
2789 * allocate. We'll allocate the return value with the same precision as
2790 * 'self'.
2791 *
2792 * Note that the caller should ensure that any argument sare removed
2793 * from the stack before calling this routine, since we might push the
2794 * self-reference onto the stack.
2795 */
setup_getp_retval(VMG_ vm_obj_id_t self,vm_val_t * retval,char ** new_ext,size_t prec)2796 int CVmObjBigNum::setup_getp_retval(VMG_ vm_obj_id_t self,
2797 vm_val_t *retval, char **new_ext,
2798 size_t prec)
2799 {
2800 /* if I'm not a number, return myself unchanged */
2801 if (is_nan(ext_))
2802 {
2803 retval->set_obj(self);
2804 return TRUE;
2805 }
2806
2807 /* push a self-reference while we're working */
2808 G_stk->push()->set_obj(self);
2809
2810 /* create a new number with the same precision as the original */
2811 retval->set_obj(create(vmg_ FALSE, prec));
2812 *new_ext = get_objid_ext(vmg_ retval->val.obj);
2813
2814 /* tell the caller to proceed */
2815 return FALSE;
2816 }
2817
2818 /*
2819 * property eval - get the fractional part
2820 */
getp_frac(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)2821 int CVmObjBigNum::getp_frac(VMG_ vm_obj_id_t self,
2822 vm_val_t *retval, uint *argc)
2823 {
2824 char *new_ext;
2825 size_t idx;
2826 int exp = get_exp(ext_);
2827 size_t prec = get_prec(ext_);
2828
2829 /* check arguments and allocate the result value */
2830 if (setup_getp_0(vmg_ self, retval, argc, &new_ext))
2831 return TRUE;
2832
2833 /* make a copy in the new object */
2834 memcpy(new_ext, ext_, calc_alloc(prec));
2835
2836 /* clear out the first n digits, where n is the exponent */
2837 for (idx = 0 ; idx < prec && (int)idx < exp ; ++idx)
2838 set_dig(new_ext, idx, 0);
2839
2840 /* normalize the result */
2841 normalize(new_ext);
2842
2843 /* remove my self-reference */
2844 G_stk->discard();
2845
2846 /* handled */
2847 return TRUE;
2848 }
2849
2850 /*
2851 * property eval - get the whole part, with no rounding
2852 */
getp_whole(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)2853 int CVmObjBigNum::getp_whole(VMG_ vm_obj_id_t self,
2854 vm_val_t *retval, uint *argc)
2855 {
2856 char *new_ext;
2857 size_t idx;
2858 int exp = get_exp(ext_);
2859 size_t prec = get_prec(ext_);
2860
2861 /* check arguments and allocate the result value */
2862 if (setup_getp_0(vmg_ self, retval, argc, &new_ext))
2863 return TRUE;
2864
2865 /* make a copy in the new object */
2866 memcpy(new_ext, ext_, calc_alloc(prec));
2867
2868 /* check what we have */
2869 if (exp <= 0)
2870 {
2871 /* it's an entirely fractional number - the result is zero */
2872 set_zero(new_ext);
2873 }
2874 else
2875 {
2876 /* clear digits after the decimal point */
2877 for (idx = (size_t)exp ; idx < prec ; ++idx)
2878 set_dig(new_ext, idx, 0);
2879
2880 /* normalize the result */
2881 normalize(new_ext);
2882 }
2883
2884 /* remove my self-reference */
2885 G_stk->discard();
2886
2887 /* handled */
2888 return TRUE;
2889 }
2890
2891 /*
2892 * property eval - round to a given number of decimal places
2893 */
getp_round_dec(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)2894 int CVmObjBigNum::getp_round_dec(VMG_ vm_obj_id_t self,
2895 vm_val_t *retval, uint *argc)
2896 {
2897 int places;
2898 char *new_ext;
2899 int post_idx;
2900 int exp = get_exp(ext_);
2901 size_t prec = get_prec(ext_);
2902 static CVmNativeCodeDesc desc(1);
2903
2904 /* check arguments */
2905 if (get_prop_check_argc(retval, argc, &desc))
2906 return TRUE;
2907
2908 /* get the number of digits for rounding */
2909 places = CVmBif::pop_int_val(vmg0_);
2910
2911 /* set up the return value */
2912 if (setup_getp_retval(vmg_ self, retval, &new_ext, prec))
2913 return TRUE;
2914
2915 /* make a copy in the new object */
2916 memcpy(new_ext, ext_, calc_alloc(prec));
2917
2918 /*
2919 * Determine if the first digit we're lopping off is actually
2920 * represented in the number or not. This digit position is the
2921 * exponent plus the number of decimal places after the decimal to
2922 * keep - if this value is at least zero and less than the
2923 * precision, it's part of the number.
2924 */
2925 post_idx = places + exp;
2926 if (post_idx < 0)
2927 {
2928 /*
2929 * the digit after the last digit to keep is actually before the
2930 * beginning of the number, so the result of the rounding is
2931 * simply zero
2932 */
2933 set_zero(new_ext);
2934 }
2935 else if (post_idx >= (int)prec)
2936 {
2937 /*
2938 * the digit after the last digit to keep is past the end of the
2939 * represented digits, so rounding has no effect at all - we'll
2940 * simply return the same number
2941 */
2942 }
2943 else
2944 {
2945 int need_to_round;
2946 size_t idx;
2947
2948 /*
2949 * the digit after the last digit is part of the number - note
2950 * it so we can tell if we need to round later
2951 */
2952 need_to_round = (get_dig(new_ext, post_idx) >= 5);
2953
2954 /* set all of the digits to be elided to zero */
2955 for (idx = (size_t)post_idx ; idx < prec ; ++idx)
2956 set_dig(new_ext, idx, 0);
2957
2958 /* if we need to round, do so now */
2959 if (need_to_round)
2960 {
2961 int carry;
2962
2963 /* increment the last digit, and apply carry as far as needed */
2964 for (carry = TRUE, idx = (size_t)post_idx ; idx != 0 ; )
2965 {
2966 /* move to the next digit */
2967 --idx;
2968
2969 /* check to see if we need to carry */
2970 if (get_dig(new_ext, idx) == 9)
2971 {
2972 /* set it to zero, then keep going to carry */
2973 set_dig(new_ext, idx, 0);
2974 }
2975 else
2976 {
2977 /* increment the digit */
2978 set_dig(new_ext, idx, get_dig(new_ext, idx) + 1);
2979
2980 /* no need to carry */
2981 carry = FALSE;
2982 break;
2983 }
2984 }
2985
2986 /* if we needed to carry, insert a leading 1 */
2987 if (carry)
2988 {
2989 /* shift the number right one place */
2990 shift_right(new_ext, 1);
2991
2992 /* adjust the exponent upwards */
2993 ++exp;
2994 set_exp(new_ext, exp);
2995
2996 /* insert the leading 1 */
2997 set_dig(new_ext, 0, 1);
2998 }
2999 }
3000
3001 /* normalize the result */
3002 normalize(new_ext);
3003 }
3004
3005 /* remove my self-reference */
3006 G_stk->discard();
3007
3008 /* handled */
3009 return TRUE;
3010 }
3011
3012 /*
3013 * property eval - get the absolute value
3014 */
getp_abs(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)3015 int CVmObjBigNum::getp_abs(VMG_ vm_obj_id_t self,
3016 vm_val_t *retval, uint *argc)
3017 {
3018 char *new_ext;
3019 size_t prec = get_prec(ext_);
3020 static CVmNativeCodeDesc desc(0);
3021
3022 /* check arguments */
3023 if (get_prop_check_argc(retval, argc, &desc))
3024 return TRUE;
3025
3026 /*
3027 * If I'm not an ordinary number or an infinity, or I'm already
3028 * non-negative, return myself unchanged. Note that we change
3029 * negative infinity to infinity, even though this might not make a
3030 * great deal of sense mathematically.
3031 */
3032 if (!get_neg(ext_)
3033 || (get_type(ext_) != VMBN_T_NUM && get_type(ext_) != VMBN_T_INF))
3034 {
3035 retval->set_obj(self);
3036 return TRUE;
3037 }
3038
3039 /* push a self-reference while we're working */
3040 G_stk->push()->set_obj(self);
3041
3042 /*
3043 * if I'm negative infinity, we don't need any precision in the new
3044 * value
3045 */
3046 if (get_type(ext_) == VMBN_T_INF)
3047 prec = 1;
3048
3049 /* create a new number with the same precision as the original */
3050 retval->set_obj(create(vmg_ FALSE, prec));
3051 new_ext = get_objid_ext(vmg_ retval->val.obj);
3052
3053 /* make a copy in the new object */
3054 memcpy(new_ext, ext_, calc_alloc(prec));
3055
3056 /* set the sign to positive */
3057 set_neg(new_ext, FALSE);
3058
3059 /* remove my self-reference */
3060 G_stk->discard();
3061
3062 /* handled */
3063 return TRUE;
3064 }
3065
3066 /*
3067 * property eval - ceiling (least integer greater than or equal to self)
3068 */
getp_ceil(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)3069 int CVmObjBigNum::getp_ceil(VMG_ vm_obj_id_t self,
3070 vm_val_t *retval, uint *argc)
3071 {
3072 char *new_ext;
3073 size_t idx;
3074 int exp = get_exp(ext_);
3075 size_t prec = get_prec(ext_);
3076 int frac_zero;
3077 static CVmNativeCodeDesc desc(0);
3078
3079 /* check arguments */
3080 if (get_prop_check_argc(retval, argc, &desc))
3081 return TRUE;
3082
3083 /* check arguments and allocate the result value */
3084 if (setup_getp_0(vmg_ self, retval, argc, &new_ext))
3085 return TRUE;
3086
3087 /* make a copy in the new object */
3088 memcpy(new_ext, ext_, calc_alloc(prec));
3089
3090 /* determine if the fractional part is non-zero */
3091 frac_zero = is_frac_zero(new_ext);
3092
3093 /* check what we have */
3094 if (exp <= 0)
3095 {
3096 /*
3097 * it's an entirely fractional number - the result is zero if
3098 * the number is negative or zero, one if the number is positive
3099 */
3100 if (get_neg(new_ext) || frac_zero)
3101 {
3102 /* -1 < x <= 0 -> ceil(x) = 0 */
3103 set_zero(new_ext);
3104 }
3105 else
3106 {
3107 /* 0 < x < 1 -> ceil(x) = 1 */
3108 copy_val(new_ext, get_one(), FALSE);
3109 }
3110 }
3111 else
3112 {
3113 /* clear digits after the decimal point */
3114 for (idx = (size_t)exp ; idx < prec ; ++idx)
3115 set_dig(new_ext, idx, 0);
3116
3117 /*
3118 * if there's a fractional part and it's positive, increment the
3119 * number
3120 */
3121 if (!frac_zero && !get_neg(new_ext))
3122 increment_abs(new_ext);
3123 }
3124
3125 /* remove my self-reference */
3126 G_stk->discard();
3127
3128 /* handled */
3129 return TRUE;
3130 }
3131
3132 /*
3133 * property eval - floor (greatest integer <= self)
3134 */
getp_floor(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)3135 int CVmObjBigNum::getp_floor(VMG_ vm_obj_id_t self,
3136 vm_val_t *retval, uint *argc)
3137 {
3138 char *new_ext;
3139 size_t idx;
3140 int exp = get_exp(ext_);
3141 size_t prec = get_prec(ext_);
3142 int frac_zero;
3143
3144 /* check arguments and allocate the result value */
3145 if (setup_getp_0(vmg_ self, retval, argc, &new_ext))
3146 return TRUE;
3147
3148 /* make a copy in the new object */
3149 memcpy(new_ext, ext_, calc_alloc(prec));
3150
3151 /* determine if the fractional part is non-zero */
3152 frac_zero = is_frac_zero(new_ext);
3153
3154 /* check what we have */
3155 if (exp <= 0)
3156 {
3157 /*
3158 * it's an entirely fractional number - the result is zero if
3159 * the number is positive or zero, -1 if the number is negative
3160 */
3161 if (!get_neg(new_ext) || frac_zero)
3162 {
3163 /* 0 <= x < 1 -> floor(x) = 0 */
3164 set_zero(new_ext);
3165 }
3166 else
3167 {
3168 /* -1 < x < 0 -> floor(x) = -1 */
3169 copy_val(new_ext, get_one(), FALSE);
3170 set_neg(new_ext, TRUE);
3171 }
3172 }
3173 else
3174 {
3175 /* clear digits after the decimal point */
3176 for (idx = (size_t)exp ; idx < prec ; ++idx)
3177 set_dig(new_ext, idx, 0);
3178
3179 /*
3180 * if there's a fractional part and the number is negative,
3181 * increment the number's absolute value
3182 */
3183 if (!frac_zero && get_neg(new_ext))
3184 increment_abs(new_ext);
3185 }
3186
3187 /* remove my self-reference */
3188 G_stk->discard();
3189
3190 /* handled */
3191 return TRUE;
3192 }
3193
3194 /*
3195 * property eval - getScale
3196 */
getp_get_scale(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)3197 int CVmObjBigNum::getp_get_scale(VMG_ vm_obj_id_t self,
3198 vm_val_t *retval, uint *argc)
3199 {
3200 static CVmNativeCodeDesc desc(0);
3201
3202 /* check arguments */
3203 if (get_prop_check_argc(retval, argc, &desc))
3204 return TRUE;
3205
3206 /* check the type */
3207 if (is_nan(ext_))
3208 {
3209 /* it's not a number - return nil for the scale */
3210 retval->set_nil();
3211 }
3212 else
3213 {
3214 /* return an integer giving the number's scale */
3215 retval->set_int(get_exp(ext_));
3216 }
3217
3218 /* handled */
3219 return TRUE;
3220 }
3221
3222 /*
3223 * property eval - scale
3224 */
getp_scale(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)3225 int CVmObjBigNum::getp_scale(VMG_ vm_obj_id_t self,
3226 vm_val_t *retval, uint *argc)
3227 {
3228 char *new_ext;
3229 size_t prec = get_prec(ext_);
3230 int scale;
3231 static CVmNativeCodeDesc desc(1);
3232
3233 /* check arguments */
3234 if (get_prop_check_argc(retval, argc, &desc))
3235 return TRUE;
3236
3237 /* get the scaling argument */
3238 scale = CVmBif::pop_int_val(vmg0_);
3239
3240 /* set up the return value */
3241 if (setup_getp_retval(vmg_ self, retval, &new_ext, prec))
3242 return TRUE;
3243
3244 /* copy the value */
3245 memcpy(new_ext, ext_, calc_alloc(prec));
3246
3247 /* adjust the exponent by the scale factor */
3248 set_exp(new_ext, get_exp(new_ext) + scale);
3249
3250 /* discard the GC protection */
3251 G_stk->discard();
3252
3253 /* handled */
3254 return TRUE;
3255 }
3256
3257 /*
3258 * property eval - negate
3259 */
getp_negate(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)3260 int CVmObjBigNum::getp_negate(VMG_ vm_obj_id_t self,
3261 vm_val_t *retval, uint *argc)
3262 {
3263 static CVmNativeCodeDesc desc(0);
3264
3265 /* check arguments */
3266 if (get_prop_check_argc(retval, argc, &desc))
3267 return TRUE;
3268
3269 /* negate the value */
3270 neg_val(vmg_ retval, self);
3271
3272 /* handled */
3273 return TRUE;
3274 }
3275
3276 /*
3277 * property eval - copy sign
3278 */
getp_copy_sign(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)3279 int CVmObjBigNum::getp_copy_sign(VMG_ vm_obj_id_t self,
3280 vm_val_t *retval, uint *argc)
3281 {
3282 vm_val_t val2;
3283 char *new_ext;
3284 const char *ext2;
3285 size_t prec = get_prec(ext_);
3286
3287 if (setup_getp_1(vmg_ self, retval, argc, &new_ext,
3288 &val2, &ext2, TRUE))
3289 return TRUE;
3290
3291 /* make a copy of my value in the new object */
3292 memcpy(new_ext, ext_, calc_alloc(prec));
3293
3294 /* set the sign from the other object */
3295 set_neg(new_ext, get_neg(ext2));
3296
3297 /*
3298 * normalize it (this is important when the value was zero to start
3299 * with, since zero is always represented without a negative sign)
3300 */
3301 normalize(new_ext);
3302
3303 /* remove the GC protection */
3304 G_stk->discard(2);
3305
3306 /* handled */
3307 return TRUE;
3308 }
3309
3310 /*
3311 * property eval - isNegative
3312 */
getp_is_neg(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)3313 int CVmObjBigNum::getp_is_neg(VMG_ vm_obj_id_t self,
3314 vm_val_t *retval, uint *argc)
3315 {
3316 static CVmNativeCodeDesc desc(0);
3317
3318 /* check arguments */
3319 if (get_prop_check_argc(retval, argc, &desc))
3320 return TRUE;
3321
3322 /* if I'm not an ordinary number or an infinity, I'm not negative */
3323 if (get_type(ext_) != VMBN_T_NUM && get_type(ext_) != VMBN_T_INF)
3324 {
3325 /* I'm not negative, so return nil */
3326 retval->set_nil();
3327 }
3328 else
3329 {
3330 /* set the return value to true if I'm negative, nil if not */
3331 retval->set_logical(get_neg(ext_));
3332 }
3333
3334 /* handled */
3335 return TRUE;
3336 }
3337
3338 /*
3339 * property eval - remainder
3340 */
getp_remainder(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)3341 int CVmObjBigNum::getp_remainder(VMG_ vm_obj_id_t self,
3342 vm_val_t *retval, uint *argc)
3343 {
3344 vm_val_t val2;
3345 const char *ext2;
3346 char *quo_ext;
3347 char *rem_ext;
3348 vm_val_t rem_val;
3349 vm_val_t quo_val;
3350 CVmObjList *lst;
3351 static CVmNativeCodeDesc desc(1);
3352
3353 /* check arguments */
3354 if (get_prop_check_argc(retval, argc, &desc))
3355 return TRUE;
3356
3357 /* pop the divisor */
3358 G_stk->pop(&val2);
3359
3360 /* convert it to BigNumber */
3361 if (!cvt_to_bignum(vmg_ self, &val2))
3362 {
3363 /* it's not a BigNumber - throw an error */
3364 err_throw(VMERR_BAD_TYPE_BIF);
3365 }
3366
3367 /* get the divisor's extension */
3368 ext2 = get_objid_ext(vmg_ val2.val.obj);
3369
3370 /* push 'self' and the other value to protect against GC */
3371 G_stk->push()->set_obj(self);
3372 G_stk->push(&val2);
3373
3374 /* create a quotient result value, and push it for safekeeping */
3375 quo_ext = compute_init_2op(vmg_ &quo_val, ext_, ext2);
3376 G_stk->push(&quo_val);
3377
3378 /* create a remainder result value, and push it for safekeeping */
3379 rem_ext = compute_init_2op(vmg_ &rem_val, ext_, ext2);
3380 G_stk->push(&rem_val);
3381
3382 /*
3383 * create a list for the return value - it will have two elements:
3384 * the quotient and the remainder
3385 */
3386 retval->set_obj(CVmObjList::create(vmg_ FALSE, 2));
3387 lst = (CVmObjList *)vm_objp(vmg_ retval->val.obj);
3388
3389 /* set the list elements */
3390 lst->cons_set_element(0, &quo_val);
3391 lst->cons_set_element(1, &rem_val);
3392
3393 /* calculate the quotient */
3394 compute_quotient_into(vmg_ quo_ext, rem_ext, ext_, ext2);
3395
3396 /* remove the GC protection */
3397 G_stk->discard(4);
3398
3399 /* handled */
3400 return TRUE;
3401 }
3402
3403 /*
3404 * property eval - sine
3405 */
getp_sin(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)3406 int CVmObjBigNum::getp_sin(VMG_ vm_obj_id_t self,
3407 vm_val_t *retval, uint *argc)
3408 {
3409 char *new_ext;
3410 size_t prec = get_prec(ext_);
3411 uint hdl1, hdl2, hdl3, hdl4, hdl5, hdl6, hdl7;
3412 char *ext1, *ext2, *ext3, *ext4, *ext5, *ext6, *ext7;
3413 int neg_result;
3414 const char *pi;
3415
3416 /* check arguments and set up the result */
3417 if (setup_getp_0(vmg_ self, retval, argc, &new_ext))
3418 return TRUE;
3419
3420 /* cache pi */
3421 pi = cache_pi(vmg_ prec + 3);
3422
3423 /*
3424 * Allocate our temporary registers. We'll use 1 and 2 to calculate
3425 * x^n - we'll start with x^(n-2) in one, and multiply by x^2 to put
3426 * the result in the other. 3 we'll use to store n!. 4 we'll use
3427 * to store the result of x^n/n!, and 5 and 6 we'll swap as the
3428 * master accumulator. 7 we'll use to store x^2.
3429 *
3430 * Allocate the temporary registers with more digits of precision
3431 * than we need in the result, to ensure that accumulated rounding
3432 * errors don't affect the result.
3433 */
3434 alloc_temp_regs(vmg_ prec + 3, 7,
3435 &ext1, &hdl1, &ext2, &hdl2, &ext3, &hdl3,
3436 &ext4, &hdl4, &ext5, &hdl5, &ext6, &hdl6, &ext7, &hdl7);
3437
3438 /* protect against errors so we definitely free the registers */
3439 err_try
3440 {
3441 char *tmp;
3442
3443 /* copy our initial value to ext1 */
3444 copy_val(ext1, ext_, FALSE);
3445
3446 /*
3447 * Note that sin(-x) = -sin(x). If x < 0, note the negative
3448 * sign and then negate the value so that we calculate sin(x)
3449 * and return the negative of the result.
3450 */
3451 neg_result = get_neg(ext1);
3452 set_neg(ext1, FALSE);
3453
3454 /* calculate 2*pi */
3455 copy_val(ext7, pi, TRUE);
3456 mul_by_long(ext7, 2);
3457
3458 /*
3459 * Because we'll use a Taylor series around 0 to calculate the
3460 * result, we want our value as close to the expansion point (0)
3461 * as possible to speed up convergence of the series.
3462 * Fortunately this is especially easy for sin() because of the
3463 * periodic nature of the function.
3464 *
3465 * First, note that sin(2*pi*i + x) = sin(x) for all integers i,
3466 * so we can reduce the argument mod 2*pi until it's in the
3467 * range 0 <= x < 2*pi (we might have to do this multiple times
3468 * if the number's scale exceeds its precision). Note that we
3469 * already made sure the number is positive.
3470 */
3471 while (compare_abs(ext1, ext7) > 0)
3472 {
3473 /* divide by 2*pi, storing the remainder in r2 */
3474 compute_quotient_into(vmg_ ext6, ext2, ext1, ext7);
3475
3476 /* swap r2 into r1 for the next round */
3477 tmp = ext1;
3478 ext1 = ext2;
3479 ext2 = tmp;
3480 }
3481
3482 /*
3483 * Next, note that sin(x+pi) = -sin(x). If x > pi, negate the
3484 * result (again if necessary) and reduce the argument by pi.
3485 * This will reduce our range to 0 <= x <= pi.
3486 */
3487 copy_val(ext7, pi, TRUE);
3488 if (compare_abs(ext1, ext7) > 0)
3489 {
3490 /* negate the result */
3491 neg_result = !neg_result;
3492
3493 /* subtract pi from the argument */
3494 compute_abs_diff_into(ext2, ext1, ext7);
3495
3496 /* swap the result into r1 */
3497 tmp = ext1;
3498 ext1 = ext2;
3499 ext2 = tmp;
3500 }
3501
3502 /*
3503 * Use the fact that sin(x + pi) = -sin(x) once again: if x >
3504 * pi/2, subtract pi from x to adjust the range to -pi/2 <= x <=
3505 * pi/2.
3506 */
3507 div_by_long(ext7, 2);
3508 if (compare_abs(ext1, ext7) > 0)
3509 {
3510 /* negate the result */
3511 neg_result = !neg_result;
3512
3513 /* subtract pi from the argument */
3514 copy_val(ext7, pi, TRUE);
3515 compute_abs_diff_into(ext2, ext1, ext7);
3516
3517 /* swap the result into r1 */
3518 tmp = ext1;
3519 ext1 = ext2;
3520 ext2 = tmp;
3521 }
3522
3523 /*
3524 * once again, reduce our range using the sign equivalence -
3525 * this will limit our range to 0 <= x <= pi/2
3526 */
3527 if (get_neg(ext1))
3528 neg_result = !neg_result;
3529 set_neg(ext1, FALSE);
3530
3531 /*
3532 * Next, observe that sin(x+pi/2) = cos(x). If x > pi/4,
3533 * calculate using the cosine series instead of the sine series.
3534 */
3535 copy_val(ext7, pi, TRUE);
3536 div_by_long(ext7, 4);
3537 if (compare_abs(ext1, ext7) > 0)
3538 {
3539 /* calculate pi/2 */
3540 copy_val(ext7, pi, TRUE);
3541 div_by_long(ext7, 2);
3542
3543 /*
3544 * subtract pi/2 - this will give us a value in the range
3545 * -pi/4 <= x <= pi/4
3546 */
3547 compute_abs_diff_into(ext2, ext1, ext7);
3548
3549 /* cos(-x) = cos(x), so we can ignore the sign */
3550 set_neg(ext1, FALSE);
3551
3552 /* swap the result into r1 */
3553 tmp = ext1;
3554 ext1 = ext2;
3555 ext2 = tmp;
3556
3557 /* calculate the cosine series */
3558 calc_cos_series(vmg_ new_ext,
3559 ext1, ext2, ext3, ext4, ext5, ext6, ext7);
3560 }
3561 else
3562 {
3563 /*
3564 * We now have a value in the range 0 <= x <= pi/4, which
3565 * will converge quickly with our Taylor series
3566 */
3567 calc_sin_series(vmg_ new_ext,
3568 ext1, ext2, ext3, ext4, ext5, ext6, ext7);
3569 }
3570
3571 /* negate the result if necessary */
3572 if (neg_result)
3573 set_neg(new_ext, !get_neg(new_ext));
3574
3575 /* normalize the result */
3576 normalize(new_ext);
3577 }
3578 err_finally
3579 {
3580 /* release our temporary registers */
3581 release_temp_regs(vmg_ 7, hdl1, hdl2, hdl3, hdl4, hdl5, hdl6, hdl7);
3582 }
3583 err_end;
3584
3585 /* remove my self-reference */
3586 G_stk->discard();
3587
3588 /* handled */
3589 return TRUE;
3590 }
3591
3592 /*
3593 * property eval - cosine. This works very much the same way as
3594 * getp_sin() - refer to the more extensive comments in that routine for
3595 * more detail on the algorithm.
3596 */
getp_cos(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)3597 int CVmObjBigNum::getp_cos(VMG_ vm_obj_id_t self,
3598 vm_val_t *retval, uint *argc)
3599 {
3600 char *new_ext;
3601 size_t prec = get_prec(ext_);
3602 uint hdl1, hdl2, hdl3, hdl4, hdl5, hdl6, hdl7;
3603 char *ext1, *ext2, *ext3, *ext4, *ext5, *ext6, *ext7;
3604 int neg_result;
3605 const char *pi;
3606
3607 /* check arguments and set up the result */
3608 if (setup_getp_0(vmg_ self, retval, argc, &new_ext))
3609 return TRUE;
3610
3611 /* cache pi to our working precision */
3612 pi = cache_pi(vmg_ prec + 3);
3613
3614 /* allocate our temporary registers, as per sin() */
3615 alloc_temp_regs(vmg_ prec + 3, 7,
3616 &ext1, &hdl1, &ext2, &hdl2, &ext3, &hdl3,
3617 &ext4, &hdl4, &ext5, &hdl5, &ext6, &hdl6, &ext7, &hdl7);
3618
3619 /* protect against errors so we definitely free the registers */
3620 err_try
3621 {
3622 char *tmp;
3623
3624 /* presume the result sign will be correct */
3625 neg_result = FALSE;
3626
3627 /* copy our initial value to ext1 */
3628 copy_val(ext1, ext_, FALSE);
3629
3630 /* note that cos(-x) = cos(x) - if x < 0, use -x */
3631 set_neg(ext1, FALSE);
3632
3633 /* reduce the argument modulo 2*pi */
3634 copy_val(ext7, pi, TRUE);
3635 mul_by_long(ext7, 2);
3636 while (compare_abs(ext1, ext7) > 0)
3637 {
3638 /* divide by 2*pi, storing the remainder in r2 */
3639 compute_quotient_into(vmg_ ext6, ext2, ext1, ext7);
3640
3641 /* swap r2 into r1 for the next round */
3642 tmp = ext1;
3643 ext1 = ext2;
3644 ext2 = tmp;
3645 }
3646
3647 /*
3648 * Next, note that cos(x+pi) = -cos(x). If x > pi, negate the
3649 * result (again if necessary) and reduce the argument by pi.
3650 * This will reduce our range to 0 <= x <= pi.
3651 */
3652 copy_val(ext7, pi, TRUE);
3653 if (compare_abs(ext1, ext7) > 0)
3654 {
3655 /* negate the result */
3656 neg_result = !neg_result;
3657
3658 /* subtract pi from the argument */
3659 compute_abs_diff_into(ext2, ext1, ext7);
3660
3661 /* swap the result into r1 */
3662 tmp = ext1;
3663 ext1 = ext2;
3664 ext2 = tmp;
3665 }
3666
3667 /*
3668 * Use the fact that cos(x + pi) = -cos(x) once again: if x >
3669 * pi/2, subtract pi from x to adjust the range to -pi/2 <= x <=
3670 * pi/2.
3671 */
3672 div_by_long(ext7, 2);
3673 if (compare_abs(ext1, ext7) > 0)
3674 {
3675 /* negate the result */
3676 neg_result = !neg_result;
3677
3678 /* subtract pi from the argument */
3679 copy_val(ext7, pi, TRUE);
3680 compute_abs_diff_into(ext2, ext1, ext7);
3681
3682 /* swap the result into r1 */
3683 tmp = ext1;
3684 ext1 = ext2;
3685 ext2 = tmp;
3686 }
3687
3688 /*
3689 * once again, reduce our range using the sign equivalence -
3690 * this will limit our range to 0 <= x <= pi/2
3691 */
3692 set_neg(ext1, FALSE);
3693
3694 /*
3695 * Next, observe that cos(x+pi/2) = -sin(x). If x > pi/4,
3696 * calculate using the sine series instead of the cosine series.
3697 */
3698 copy_val(ext7, pi, TRUE);
3699 div_by_long(ext7, 4);
3700 if (compare_abs(ext1, ext7) > 0)
3701 {
3702 /* negate the result */
3703 neg_result = !neg_result;
3704
3705 /* calculate pi/2 */
3706 copy_val(ext7, pi, TRUE);
3707 div_by_long(ext7, 2);
3708
3709 /*
3710 * subtract pi/2 - this will give us a value in the range
3711 * -pi/4 <= x <= pi/4
3712 */
3713 compute_abs_diff_into(ext2, ext1, ext7);
3714
3715 /* sin(-x) = -sin(x) */
3716 if (get_neg(ext1))
3717 neg_result = !neg_result;
3718 set_neg(ext1, FALSE);
3719
3720 /* swap the result into r1 */
3721 tmp = ext1;
3722 ext1 = ext2;
3723 ext2 = tmp;
3724
3725 /* calculate the sine series */
3726 calc_sin_series(vmg_ new_ext,
3727 ext1, ext2, ext3, ext4, ext5, ext6, ext7);
3728 }
3729 else
3730 {
3731 /*
3732 * We now have a value in the range 0 <= x <= pi/4, which
3733 * will converge quickly with our Taylor series
3734 */
3735 calc_cos_series(vmg_ new_ext,
3736 ext1, ext2, ext3, ext4, ext5, ext6, ext7);
3737 }
3738
3739 /* negate the result if necessary */
3740 if (neg_result)
3741 set_neg(new_ext, !get_neg(new_ext));
3742
3743 /* normalize the result */
3744 normalize(new_ext);
3745 }
3746 err_finally
3747 {
3748 /* release our temporary registers */
3749 release_temp_regs(vmg_ 7, hdl1, hdl2, hdl3, hdl4, hdl5, hdl6, hdl7);
3750 }
3751 err_end;
3752
3753 /* remove my self-reference */
3754 G_stk->discard();
3755
3756 /* handled */
3757 return TRUE;
3758 }
3759
3760 /*
3761 * property eval - tangent. We calculate the sine and cosine, then
3762 * compute the quotient to determine the tangent. This routine works
3763 * very much like the sin() and cos() routines; refer to getp_sin() for
3764 * more thorough documentation.
3765 */
getp_tan(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)3766 int CVmObjBigNum::getp_tan(VMG_ vm_obj_id_t self,
3767 vm_val_t *retval, uint *argc)
3768 {
3769 char *new_ext;
3770 size_t prec = get_prec(ext_);
3771 uint hdl1, hdl2, hdl3, hdl4, hdl5, hdl6, hdl7, hdl8, hdl9;
3772 char *ext1, *ext2, *ext3, *ext4, *ext5, *ext6, *ext7, *ext8, *ext9;
3773 int neg_result;
3774 int invert_result;
3775 const char *pi;
3776
3777 /* check arguments and set up the result */
3778 if (setup_getp_0(vmg_ self, retval, argc, &new_ext))
3779 return TRUE;
3780
3781 /* cache pi to the working precision */
3782 pi = cache_pi(vmg_ prec + 3);
3783
3784 /*
3785 * Allocate our temporary registers for sin() and cos()
3786 * calculations, plus two extra: one to hold the sine while we're
3787 * calculating the cosine, and the other to hold the cosine result.
3788 *
3789 * (We could make do with fewer registers by copying values around,
3790 * but if the numbers are of very high precision it's much cheaper
3791 * to allocate more registers, since the registers come out of the
3792 * register cache and probably won't require any actual memory
3793 * allocation.)
3794 */
3795 alloc_temp_regs(vmg_ prec + 3, 9,
3796 &ext1, &hdl1, &ext2, &hdl2, &ext3, &hdl3,
3797 &ext4, &hdl4, &ext5, &hdl5, &ext6, &hdl6,
3798 &ext7, &hdl7, &ext8, &hdl8, &ext9, &hdl9);
3799
3800 /* protect against errors so we definitely free the registers */
3801 err_try
3802 {
3803 char *tmp;
3804
3805 /* presume we won't have to invert our result */
3806 invert_result = FALSE;
3807
3808 /* copy our initial value to ext1 */
3809 copy_val(ext1, ext_, FALSE);
3810
3811 /* note that tan(-x) = -tan(x) - if x < 0, use -x */
3812 neg_result = get_neg(ext1);
3813 set_neg(ext1, FALSE);
3814
3815 /* reduce the argument modulo 2*pi */
3816 copy_val(ext7, pi, TRUE);
3817 mul_by_long(ext7, 2);
3818 while (compare_abs(ext1, ext7) > 0)
3819 {
3820 /* divide by 2*pi, storing the remainder in r2 */
3821 compute_quotient_into(vmg_ ext6, ext2, ext1, ext7);
3822
3823 /* swap r2 into r1 for the next round */
3824 tmp = ext1;
3825 ext1 = ext2;
3826 ext2 = tmp;
3827 }
3828
3829 /*
3830 * Next, note that tan(x+pi) = tan(x). So, if x > pi, the
3831 * argument by pi. This will reduce our range to 0 <= x <= pi.
3832 */
3833 copy_val(ext7, pi, TRUE);
3834 if (compare_abs(ext1, ext7) > 0)
3835 {
3836 /* subtract pi from the argument */
3837 compute_abs_diff_into(ext2, ext1, ext7);
3838
3839 /* swap the result into r1 */
3840 tmp = ext1;
3841 ext1 = ext2;
3842 ext2 = tmp;
3843 }
3844
3845 /*
3846 * Use the fact that tan(x + pi) = tan(x) once again: if x >
3847 * pi/2, subtract pi from x to adjust the range to -pi/2 <= x <=
3848 * pi/2.
3849 */
3850 div_by_long(ext7, 2);
3851 if (compare_abs(ext1, ext7) > 0)
3852 {
3853 /* subtract pi from the argument */
3854 copy_val(ext7, pi, TRUE);
3855 compute_abs_diff_into(ext2, ext1, ext7);
3856
3857 /* swap the result into r1 */
3858 tmp = ext1;
3859 ext1 = ext2;
3860 ext2 = tmp;
3861 }
3862
3863 /*
3864 * once again, reduce our range using the sign equivalence -
3865 * this will limit our range to 0 <= x <= pi/2
3866 */
3867 if (get_neg(ext1))
3868 neg_result = !neg_result;
3869 set_neg(ext1, FALSE);
3870
3871 /*
3872 * Next, observe that tan(x+pi/2) = 1/tan(x). If x > pi/4,
3873 * invert the result.
3874 */
3875 copy_val(ext7, pi, TRUE);
3876 div_by_long(ext7, 4);
3877 if (compare_abs(ext1, ext7) > 0)
3878 {
3879 /* calculate pi/2 */
3880 copy_val(ext7, pi, TRUE);
3881 div_by_long(ext7, 2);
3882
3883 /* subtract pi/2 to get into range -pi/4 <= x <= pi/4 */
3884 compute_abs_diff_into(ext2, ext1, ext7);
3885
3886 /* sin(-x) = -sin(x) */
3887 if (get_neg(ext1))
3888 neg_result = !neg_result;
3889 set_neg(ext1, FALSE);
3890
3891 /* swap the result into r1 */
3892 tmp = ext1;
3893 ext1 = ext2;
3894 ext2 = tmp;
3895
3896 /* note that we must invert the result */
3897 invert_result = TRUE;
3898 }
3899
3900 /*
3901 * make a copy of our argument value in ext9 for safekeeping
3902 * while we're calculating the sine
3903 */
3904 copy_val(ext9, ext1, FALSE);
3905
3906 /*
3907 * We now have a value in the range 0 <= x <= pi/4, which will
3908 * converge quickly with our Taylor series for sine and cosine.
3909 * This will also ensure that both sin and cos are non-negative,
3910 * so the sign we calculated for the tangent is all that matters
3911 * for the sign.
3912 *
3913 * First, Calculate the sine and store the result in r8.
3914 */
3915 calc_sin_series(vmg_ ext8,
3916 ext1, ext2, ext3, ext4, ext5, ext6, ext7);
3917
3918 /*
3919 * Calculate the cosine and store the result in r1. Note that
3920 * we saved the argument value in ext9 while we were working on
3921 * the sine, so we can now use that value as the argument here.
3922 * ext1 was trashed by the sine calculation, so just use it as
3923 * the output register here.
3924 */
3925 calc_cos_series(vmg_ ext1,
3926 ext9, ext2, ext3, ext4, ext5, ext6, ext7);
3927
3928 /* calculate the quotient sin/cos, or cos/sin if inverted */
3929 if (invert_result)
3930 compute_quotient_into(vmg_ new_ext, 0, ext1, ext8);
3931 else
3932 compute_quotient_into(vmg_ new_ext, 0, ext8, ext1);
3933
3934 /* negate the result if necessary */
3935 set_neg(new_ext, neg_result);
3936
3937 /* normalize the result */
3938 normalize(new_ext);
3939 }
3940 err_finally
3941 {
3942 /* release our temporary registers */
3943 release_temp_regs(vmg_ 9, hdl1, hdl2, hdl3, hdl4,
3944 hdl5, hdl6, hdl7, hdl8, hdl9);
3945 }
3946 err_end;
3947
3948 /* remove my self-reference */
3949 G_stk->discard();
3950
3951 /* handled */
3952 return TRUE;
3953 }
3954
3955 /*
3956 * property evaluator - convert degrees to radians
3957 */
getp_deg2rad(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)3958 int CVmObjBigNum::getp_deg2rad(VMG_ vm_obj_id_t self,
3959 vm_val_t *retval, uint *argc)
3960 {
3961 char *new_ext;
3962 size_t prec = get_prec(ext_);
3963 uint hdl1;
3964 char *ext1;
3965 const char *pi;
3966
3967 /* check arguments and set up the result */
3968 if (setup_getp_0(vmg_ self, retval, argc, &new_ext))
3969 return TRUE;
3970
3971 /* cache pi to our required precision */
3972 pi = cache_pi(vmg_ prec + 2);
3973
3974 /*
3975 * allocate a temporary register for pi/180 - give it some extra
3976 * precision
3977 */
3978 alloc_temp_regs(vmg_ prec + 2, 1, &ext1, &hdl1);
3979
3980 /* get pi to our precision */
3981 copy_val(ext1, pi, TRUE);
3982
3983 /* divide pi by 180 */
3984 div_by_long(ext1, 180);
3985
3986 /* go back to our working precision, rounding if necessary */
3987 set_prec(ext1, prec);
3988 if (get_dig(ext1, prec) >= 5)
3989 round_up_abs(ext1);
3990
3991 /* multiply our value by pi/180 */
3992 compute_prod_into(new_ext, ext_, ext1);
3993
3994 /* release our temporary registers */
3995 release_temp_regs(vmg_ 1, hdl1);
3996
3997 /* remove my self-reference */
3998 G_stk->discard();
3999
4000 /* handled */
4001 return TRUE;
4002 }
4003
4004 /*
4005 * property evaluator - convert radians to degrees
4006 */
getp_rad2deg(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)4007 int CVmObjBigNum::getp_rad2deg(VMG_ vm_obj_id_t self,
4008 vm_val_t *retval, uint *argc)
4009 {
4010 char *new_ext;
4011 size_t prec = get_prec(ext_);
4012 uint hdl1;
4013 char *ext1;
4014 const char *pi;
4015
4016 /* check arguments and set up the result */
4017 if (setup_getp_0(vmg_ self, retval, argc, &new_ext))
4018 return TRUE;
4019
4020 /* cache pi to our working precision */
4021 pi = cache_pi(vmg_ prec + 2);
4022
4023 /* allocate a temporary register for pi/180 */
4024 alloc_temp_regs(vmg_ prec + 2, 1, &ext1, &hdl1);
4025
4026 /* catch errors so we can be sure to free registers */
4027 err_try
4028 {
4029 /* get pi to our precision */
4030 copy_val(ext1, pi, TRUE);
4031
4032 /* divide pi by 180 */
4033 div_by_long(ext1, 180);
4034
4035 /* go back to our working precision, rounding if necessary */
4036 set_prec(ext1, prec);
4037 if (get_dig(ext1, prec) >= 5)
4038 round_up_abs(ext1);
4039
4040 /* divide by pi/180 */
4041 compute_quotient_into(vmg_ new_ext, 0, ext_, ext1);
4042 }
4043 err_finally
4044 {
4045 /* release our temporary registers */
4046 release_temp_regs(vmg_ 1, hdl1);
4047 }
4048 err_end;
4049
4050 /* remove my self-reference */
4051 G_stk->discard();
4052
4053 /* handled */
4054 return TRUE;
4055 }
4056
4057 /*
4058 * property evaluator - arcsine
4059 */
getp_asin(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)4060 int CVmObjBigNum::getp_asin(VMG_ vm_obj_id_t self,
4061 vm_val_t *retval, uint *argc)
4062 {
4063 /* calculate and return the arcsine */
4064 return calc_asincos(vmg_ self, retval, argc, FALSE);
4065 }
4066
4067 /*
4068 * property evaluator - arccosine
4069 */
getp_acos(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)4070 int CVmObjBigNum::getp_acos(VMG_ vm_obj_id_t self,
4071 vm_val_t *retval, uint *argc)
4072 {
4073 /* calculate and return the arcsine */
4074 return calc_asincos(vmg_ self, retval, argc, TRUE);
4075 }
4076
4077 /*
4078 * Common property evaluator routine - arcsine and arccosine. arcsin
4079 * and arccos are related by arccos(x) = pi/2 - arcsin(x). So, to
4080 * calculate an arccos, we can simply calculate the arcsin, then
4081 * subtract the result from pi/2.
4082 */
calc_asincos(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc,int is_acos)4083 int CVmObjBigNum::calc_asincos(VMG_ vm_obj_id_t self,
4084 vm_val_t *retval, uint *argc, int is_acos)
4085 {
4086 char *new_ext;
4087
4088 /* check arguments and set up the result */
4089 if (setup_getp_0(vmg_ self, retval, argc, &new_ext))
4090 return TRUE;
4091
4092 /* calculate the arcsin/arccos into the result */
4093 calc_asincos_into(vmg_ new_ext, ext_, is_acos);
4094
4095 /* remove my self-reference */
4096 G_stk->discard();
4097
4098 /* handled */
4099 return TRUE;
4100 }
4101
4102 /*
4103 * Calculate the arcsine or arccosine into the given result variable
4104 */
calc_asincos_into(VMG_ char * dst,const char * src,int is_acos)4105 void CVmObjBigNum::calc_asincos_into(VMG_ char *dst, const char *src,
4106 int is_acos)
4107 {
4108 size_t prec = get_prec(dst);
4109 uint hdl1, hdl2, hdl3, hdl4, hdl5;
4110 char *ext1, *ext2, *ext3, *ext4, *ext5;
4111 const char *pi;
4112
4113 /* cache pi to our working precision */
4114 pi = cache_pi(vmg_ prec + 3);
4115
4116 /*
4117 * allocate our temporary registers - use some extra precision over
4118 * what we need for the result, to reduce the effect of accumulated
4119 * rounding error
4120 */
4121 alloc_temp_regs(vmg_ prec + 3, 5,
4122 &ext1, &hdl1, &ext2, &hdl2, &ext3, &hdl3,
4123 &ext4, &hdl4, &ext5, &hdl5);
4124
4125 /* catch errors so we release our temp registers */
4126 err_try
4127 {
4128 char *tmp;
4129 static const unsigned char one_over_sqrt2[] =
4130 {
4131 /* precision = 10, scale = 0, flags = 0 */
4132 10, 0, 0, 0, 0,
4133 0x70, 0x71, 0x06, 0x78, 0x14
4134 };
4135 int use_sqrt;
4136 int sqrt_neg;
4137
4138 /* get the initial value of x into our accumulator, r1 */
4139 copy_val(ext1, src, FALSE);
4140
4141 /*
4142 * check to see if the absolute value of the argument is greater
4143 * than 1 - if it is, it's not valid
4144 */
4145 copy_val(ext2, get_one(), FALSE);
4146 if (compare_abs(ext1, ext2) > 0)
4147 err_throw(VMERR_OUT_OF_RANGE);
4148
4149 /* presume we won't need to use the sqrt(1-x^2) method */
4150 use_sqrt = FALSE;
4151
4152 /*
4153 * Check to see if the absolute value of the argument is greater
4154 * than 1/sqrt(2). If it is, the series expansion converges too
4155 * slowly (as slowly as never, if the value is exactly 1). To
4156 * speed things up, in these cases calculate pi/2 -
4157 * asin(sqrt(1-x^2)) instead, which is equivalent but gives us a
4158 * smaller asin() argument for quicker series convergence.
4159 *
4160 * We don't need to compare to 1/sqrt(2) in great precision;
4161 * just use a few digits.
4162 */
4163 copy_val(ext2, (const char *)one_over_sqrt2, TRUE);
4164 if (compare_abs(ext1, ext2) > 0)
4165 {
4166 /* flag that we're using the sqrt method */
4167 use_sqrt = TRUE;
4168
4169 /* note the sign - we'll need to apply this to the result */
4170 sqrt_neg = get_neg(ext1);
4171
4172 /* compute x^2 into r2 */
4173 compute_prod_into(ext2, ext1, ext1);
4174
4175 /* subtract r2 from 1 (by adding -r2 to 1), storing in r4 */
4176 copy_val(ext3, get_one(), FALSE);
4177 make_negative(ext2);
4178 compute_sum_into(ext4, ext3, ext2);
4179
4180 /* compute sqrt(1-x^2) (which is sqrt(r4)) into r1 */
4181 compute_sqrt_into(vmg_ ext1, ext4);
4182 }
4183
4184 /* compute the arcsine */
4185 ext1 = calc_asin_series(ext1, ext2, ext3, ext4, ext5);
4186
4187 /* if we're using the sqrt method, finish the sqrt calculation */
4188 if (use_sqrt)
4189 {
4190 /* calculate pi/2 */
4191 copy_val(ext2, pi, TRUE);
4192 div_by_long(ext2, 2);
4193
4194 /* compute pi/2 - r1 by negating r1 and adding it */
4195 negate(ext1);
4196 compute_sum_into(ext3, ext2, ext1);
4197
4198 /* negate the result if the original value was negative */
4199 if (sqrt_neg)
4200 negate(ext3);
4201
4202 /* swap the result back into r1 */
4203 tmp = ext1;
4204 ext1 = ext3;
4205 ext3 = tmp;
4206 }
4207
4208 /*
4209 * We now have the arcsine in r1. If we actually wanted the
4210 * arccosine, subtract the arcsine from pi/2 to yield the
4211 * arccosine.
4212 */
4213 if (is_acos)
4214 {
4215 /* get pi/2 into r2 */
4216 copy_val(ext2, pi, TRUE);
4217 div_by_long(ext2, 2);
4218
4219 /* negate r1 to get -asin */
4220 negate(ext1);
4221
4222 /* add -asin to r2 to yield the arccosine in r3 */
4223 compute_sum_into(ext3, ext2, ext1);
4224
4225 /* swap the result back into ext1 */
4226 tmp = ext3;
4227 ext3 = ext1;
4228 ext1 = tmp;
4229 }
4230
4231 /* store the result, rounding if necessary */
4232 copy_val(dst, ext1, TRUE);
4233 }
4234 err_finally
4235 {
4236 /* release our temporary registers */
4237 release_temp_regs(vmg_ 5, hdl1, hdl2, hdl3, hdl4, hdl5);
4238 }
4239 err_end;
4240 }
4241
4242 /*
4243 * Calculate the asin series expansion. This should only be called with
4244 * argument values with absolute value less than 1/sqrt(2), because the
4245 * series converges very slowly for larger values. For operands above
4246 * 1/sqrt(2), the caller should instead compute the equivalent value
4247 *
4248 * +/- (pi/2 - asin(sqrt(1-x^2))).
4249 *
4250 * (+ if x > 0, - if x < 0).
4251 *
4252 * The argument value is in ext1, and we return a pointer to the
4253 * register that contains the result (which will be one of the input
4254 * registers). We use all of the input registers as scratchpads, so
4255 * their values are not retained.
4256 */
calc_asin_series(char * ext1,char * ext2,char * ext3,char * ext4,char * ext5)4257 char *CVmObjBigNum::calc_asin_series(char *ext1, char *ext2,
4258 char *ext3, char *ext4, char *ext5)
4259 {
4260 ulong n;
4261
4262 /* get the current power of x (1) into the x power register, r2 */
4263 copy_val(ext2, ext1, FALSE);
4264
4265 /*
4266 * compute x^2 into r3 - we'll multiply the previous power by this
4267 * to get the next power (we need x^1, x^3, x^5, etc)
4268 */
4269 compute_prod_into(ext3, ext1, ext1);
4270
4271 /* start at the first term */
4272 n = 1;
4273
4274 /* keep going until we have enough precision in the result */
4275 for (;;)
4276 {
4277 ulong i;
4278 char *tmp;
4279
4280 /* move on to the next term */
4281 n += 2;
4282
4283 /*
4284 * compute the next weirdness factor into r4:
4285 *
4286 * 1*3*5*...*(n-2)
4287 *. -----------------
4288 *. 2*4*6*...*(n-1)*n
4289 */
4290
4291 /* start out with 1 in r4 */
4292 copy_val(ext4, get_one(), FALSE);
4293
4294 /* multiply by odd numbers up to but not including 'n' */
4295 for (i = 3 ; i < n ; i += 2)
4296 mul_by_long(ext4, i);
4297
4298 /* divide by even numbers up to and including n-1 */
4299 for (i = 2 ; i < n ; i += 2)
4300 div_by_long(ext4, i);
4301
4302 /* divide by n */
4303 div_by_long(ext4, n);
4304
4305 /*
4306 * compute the next power of x - multiply our current power of x
4307 * (r2) by x^2 (r3)
4308 */
4309 compute_prod_into(ext5, ext2, ext3);
4310
4311 /* swap r5 into r2 - this new power of x is now current */
4312 tmp = ext5;
4313 ext5 = ext2;
4314 ext2 = tmp;
4315
4316 /*
4317 * multiply the current x power by the current weirdness factor
4318 * - this will yield the current term into r5
4319 */
4320 compute_prod_into(ext5, ext2, ext4);
4321
4322 /*
4323 * if this value is too small to show up in our accumulator,
4324 * we're done
4325 */
4326 if (is_zero(ext5)
4327 || get_exp(ext1) - get_exp(ext5) > (int)get_prec(ext1))
4328 break;
4329
4330 /*
4331 * we can trash the weird factor now - use it as a scratch pad
4332 * for adding the accumulator so far (r1) to this term
4333 */
4334 compute_sum_into(ext4, ext1, ext5);
4335
4336 /* swap the result into r1, since it's the new accumulator */
4337 tmp = ext4;
4338 ext4 = ext1;
4339 ext1 = tmp;
4340 }
4341
4342 /* return the accumulator register */
4343 return ext1;
4344 }
4345
4346 /*
4347 * property evaluator - arctangent
4348 */
getp_atan(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)4349 int CVmObjBigNum::getp_atan(VMG_ vm_obj_id_t self,
4350 vm_val_t *retval, uint *argc)
4351 {
4352 char *new_ext;
4353 size_t prec = get_prec(ext_);
4354 uint hdl1, hdl2, hdl3, hdl4, hdl5;
4355 char *ext1, *ext2, *ext3, *ext4, *ext5;
4356 const char *pi;
4357
4358 /* check arguments and set up the result */
4359 if (setup_getp_0(vmg_ self, retval, argc, &new_ext))
4360 return TRUE;
4361
4362 /* cache pi to our working precision */
4363 pi = cache_pi(vmg_ prec + 3);
4364
4365 /* allocate some temporary registers */
4366 alloc_temp_regs(vmg_ prec + 3, 5,
4367 &ext1, &hdl1, &ext2, &hdl2, &ext3, &hdl3,
4368 &ext4, &hdl4, &ext5, &hdl5);
4369
4370 /* catch errors so we can be sure to free registers */
4371 err_try
4372 {
4373 /*
4374 * If x has absolute value close to 1, either above or below,
4375 * our two series don't converge very rapidly. Happily, we can
4376 * fall back on an alternative in these case by using the
4377 * relation
4378 *
4379 * arctan(x) = (+/-)arccos(1/sqrt(x^2+1))
4380 *
4381 * (The sign of the result is the same as the sign of x.) Since
4382 * we already have routines for arccosine and square root, this
4383 * calculation is easy.
4384 *
4385 * If x doesn't have absolute value close to 1, use the
4386 * appropriate series, since they converge rapidly.
4387 */
4388 if (get_exp(ext_) < -1 || get_exp(ext_) > 2)
4389 {
4390 int term_neg;
4391 ulong n;
4392
4393 /*
4394 * exponent less than -1 -> the number has a small absolute
4395 * value - use the small-x series expansion:
4396 *
4397 * x - x^3/3 + x^5/5 - x^7/7 ...
4398 *
4399 * OR
4400 *
4401 * exponent greater than 2 -> the number has a large
4402 * absolute value, so the large-x series expansion should
4403 * converge quickly:
4404 *
4405 * +/- pi/2 - 1/x + 1/3x^3 - 1/5x^5 ...
4406 *
4407 * (the sign of the leading pi/2 term is positive if x is
4408 * positive, negative if x is negative)
4409 *
4410 * We can commonize these expressions by defining x' = x for
4411 * small x and x' = 1/x for large x, defining C as 0 for
4412 * small x and +/-pi/2 for large X, defining N as +1 for
4413 * small x and -1 for large x, and rewriting the series as:
4414 *
4415 * C + Nx' - Nx'^3/3 + Nx'^5/5 + ...
4416 */
4417
4418 /* check for large or small value */
4419 if (get_exp(ext_) < 0)
4420 {
4421 /* small number - start with zero in the accumulator (r1) */
4422 set_zero(ext1);
4423
4424 /* get the current power of x' into r2 - this is simply x */
4425 copy_val(ext2, ext_, FALSE);
4426
4427 /* the first term (x) is positive */
4428 term_neg = FALSE;
4429 }
4430 else
4431 {
4432 /* large number - start with pi/2 in the accumulator (r1) */
4433 copy_val(ext1, pi, TRUE);
4434 div_by_long(ext1, 2);
4435
4436 /* if x is negative, make that -pi/2 */
4437 set_neg(ext1, get_neg(ext_));
4438
4439 /* get 1/x into r2 - this is our x' term value */
4440 compute_quotient_into(vmg_ ext2, 0, get_one(), ext_);
4441
4442 /* the first term (1/x) is negative */
4443 term_neg = TRUE;
4444 }
4445
4446 /* start at the first term */
4447 n = 1;
4448
4449 /*
4450 * get x'^2 into r3 - we'll use this to calculate each
4451 * subsequent term (we need x', x'^3, x'^5, etc)
4452 */
4453 compute_prod_into(ext3, ext2, ext2);
4454
4455 /* iterate until we reach the desired precision */
4456 for (;;)
4457 {
4458 char *tmp;
4459
4460 /* copy the current power term from r2 into r4 */
4461 copy_val(ext4, ext2, FALSE);
4462
4463 /* divide by the current term constant */
4464 div_by_long(ext4, n);
4465
4466 /* negate this term if necessary */
4467 if (term_neg)
4468 set_neg(ext4, !get_neg(ext4));
4469
4470 /*
4471 * if we're under the radar on precision, stop looping
4472 * (don't stop on the first term, since the accumulator
4473 * hasn't been fully primed yet)
4474 */
4475 if (n != 1
4476 && (is_zero(ext4)
4477 || (get_exp(ext1) - get_exp(ext4)
4478 > (int)get_prec(ext1))))
4479 break;
4480
4481 /* compute the sum of the accumulator and this term */
4482 compute_sum_into(ext5, ext1, ext4);
4483
4484 /* swap the result back into the accumulator */
4485 tmp = ext1;
4486 ext1 = ext5;
4487 ext5 = tmp;
4488
4489 /*
4490 * move on to the next term - advance n by 2 and swap
4491 * the term sign
4492 */
4493 n += 2;
4494 term_neg = !term_neg;
4495
4496 /*
4497 * compute the next power term - multiply r2 (the latest
4498 * x' power) by r3 (x'^2)
4499 */
4500 compute_prod_into(ext4, ext2, ext3);
4501
4502 /* swap r4 back into r2 - this is now the latest power */
4503 tmp = ext2;
4504 ext2 = ext4;
4505 ext4 = tmp;
4506 }
4507
4508 /* store the result, rounding as needed */
4509 copy_val(new_ext, ext1, TRUE);
4510 }
4511 else
4512 {
4513 /*
4514 * We have a value of x from .01 to 10 - in this range, the
4515 * rewrite using arccosine will give us excellent precision.
4516 */
4517
4518 /* calculate x^2 into r1 */
4519 compute_prod_into(ext1, ext_, ext_);
4520
4521 /* add one (x^2 has to be positive, so don't worry about sign) */
4522 increment_abs(ext1);
4523
4524 /* take the square root and put the result in r2 */
4525 compute_sqrt_into(vmg_ ext2, ext1);
4526
4527 /* divide that into 1, and put the result back in r1 */
4528 compute_quotient_into(vmg_ ext1, 0, get_one(), ext2);
4529
4530 /*
4531 * Compute the arccosine of this value - the result is the
4532 * arctangent, so we can store it directly in the output
4533 * register.
4534 */
4535 calc_asincos_into(vmg_ new_ext, ext1, TRUE);
4536
4537 /* if the input was negative, invert the sign of the result */
4538 if (get_neg(ext_))
4539 negate(new_ext);
4540 }
4541 }
4542 err_finally
4543 {
4544 /* release our temporary registers */
4545 release_temp_regs(vmg_ 5, hdl1, hdl2, hdl3, hdl4, hdl5);
4546 }
4547 err_end;
4548
4549 /* discard the GC protection */
4550 G_stk->discard();
4551
4552 /* handled */
4553 return TRUE;
4554 }
4555
4556 /*
4557 * property evaluator - square root
4558 */
getp_sqrt(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)4559 int CVmObjBigNum::getp_sqrt(VMG_ vm_obj_id_t self,
4560 vm_val_t *retval, uint *argc)
4561 {
4562 char *new_ext;
4563
4564 /* check arguments and set up the result */
4565 if (setup_getp_0(vmg_ self, retval, argc, &new_ext))
4566 return TRUE;
4567
4568 /* compute the square root into the result */
4569 compute_sqrt_into(vmg_ new_ext, ext_);
4570
4571 /* discard the GC protection */
4572 G_stk->discard();
4573
4574 /* handled */
4575 return TRUE;
4576 }
4577
4578 /*
4579 * property evaluator - natural logarithm
4580 */
getp_ln(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)4581 int CVmObjBigNum::getp_ln(VMG_ vm_obj_id_t self,
4582 vm_val_t *retval, uint *argc)
4583 {
4584 char *new_ext;
4585
4586 /* zero or negative values are not allowed */
4587 if (is_zero(ext_) || get_neg(ext_))
4588 err_throw(VMERR_OUT_OF_RANGE);
4589
4590 /* check arguments and set up the result */
4591 if (setup_getp_0(vmg_ self, retval, argc, &new_ext))
4592 return TRUE;
4593
4594 /* compute the natural logarithm */
4595 compute_ln_into(vmg_ new_ext, ext_);
4596
4597 /* discard the GC protection */
4598 G_stk->discard();
4599
4600 /* handled */
4601 return TRUE;
4602 }
4603
4604 /*
4605 * property evaluator - exp (exponentiate e, the base of the natural
4606 * logarithm, to the power of this number)
4607 */
getp_exp(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)4608 int CVmObjBigNum::getp_exp(VMG_ vm_obj_id_t self,
4609 vm_val_t *retval, uint *argc)
4610 {
4611 char *new_ext;
4612
4613 /* check arguments and set up the result */
4614 if (setup_getp_0(vmg_ self, retval, argc, &new_ext))
4615 return TRUE;
4616
4617 /* compute exp(x) */
4618 compute_exp_into(vmg_ new_ext, ext_);
4619
4620 /* discard the GC protection */
4621 G_stk->discard();
4622
4623 /* handled */
4624 return TRUE;
4625 }
4626
4627 /*
4628 * property evaluator - log10 (base-10 logarithm)
4629 */
getp_log10(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)4630 int CVmObjBigNum::getp_log10(VMG_ vm_obj_id_t self,
4631 vm_val_t *retval, uint *argc)
4632 {
4633 char *new_ext;
4634 size_t prec = get_prec(ext_);
4635 uint hdl1, hdl2, hdl3;
4636 char *ext1, *ext2, *ext3;
4637 const char *ln10;
4638
4639 /* check arguments and set up the result */
4640 if (setup_getp_0(vmg_ self, retval, argc, &new_ext))
4641 return TRUE;
4642
4643 /* cache ln(10) to the required precision */
4644 ln10 = cache_ln10(vmg_ prec + 3);
4645
4646 /* allocate some temporary registers */
4647 alloc_temp_regs(vmg_ prec + 3, 3,
4648 &ext1, &hdl1, &ext2, &hdl2, &ext3, &hdl3);
4649
4650 /* catch errors so we can be sure to free registers */
4651 err_try
4652 {
4653 /* compute the natural logarithm of the number */
4654 compute_ln_into(vmg_ ext1, ext_);
4655
4656 /* get ln(10), properly rounded */
4657 copy_val(ext2, ln10, TRUE);
4658
4659 /* compute ln(x)/ln(10) to yield log10(x) */
4660 compute_quotient_into(vmg_ ext3, 0, ext1, ext2);
4661
4662 /* store the result, rounding as needed */
4663 copy_val(new_ext, ext3, TRUE);
4664 }
4665 err_finally
4666 {
4667 /* release our temporary registers */
4668 release_temp_regs(vmg_ 3, hdl1, hdl2, hdl3);
4669 }
4670 err_end;
4671
4672 /* discard the GC protection */
4673 G_stk->discard();
4674
4675 /* handled */
4676 return TRUE;
4677 }
4678
4679 /*
4680 * property evaluator - pow(y) - calculate x^y
4681 */
getp_pow(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)4682 int CVmObjBigNum::getp_pow(VMG_ vm_obj_id_t self,
4683 vm_val_t *retval, uint *argc)
4684 {
4685 vm_val_t val2;
4686 const char *val2_ext;
4687 char *new_ext;
4688 size_t prec;
4689 uint hdl1, hdl2;
4690 char *ext1, *ext2;
4691
4692 /* check arguments and allocate the result value */
4693 if (setup_getp_1(vmg_ self, retval, argc, &new_ext,
4694 &val2, &val2_ext, FALSE))
4695 return TRUE;
4696
4697 /* use the precision of the result */
4698 prec = get_prec(new_ext);
4699
4700 /*
4701 * Check for a special case: if the number we're exponentiating is
4702 * zero, the result is 1 for any exponent. This is a special case
4703 * because our general-purpose calculation requires that we take the
4704 * log of the base, but we can't take the log of zero.
4705 */
4706 if (is_zero(ext_))
4707 {
4708 /* 0^0 is undefined */
4709 if (is_zero(val2_ext))
4710 err_throw(VMERR_OUT_OF_RANGE);
4711
4712 /* set the result to one, and we're done */
4713 copy_val(new_ext, get_one(), FALSE);
4714 goto done;
4715 }
4716
4717 /* allocate some temporary registers */
4718 alloc_temp_regs(vmg_ prec + 3, 2,
4719 &ext1, &hdl1, &ext2, &hdl2);
4720
4721 /* catch errors so we can be sure to free registers */
4722 err_try
4723 {
4724 int result_neg;
4725
4726 /*
4727 * If a = e^b, then b = ln(a). This means that x^y = e^ln(x^y)
4728 * = e^(y * ln(x)). So, we can compute the result in terms of
4729 * natural logarithm and exponentiation of 'e', for which we
4730 * have primitives we can call.
4731 */
4732
4733 /*
4734 * If x is negative, we can only exponentiate the value to
4735 * integer powers. In this case, we can substitute x' = -x
4736 * (hence x' will be positive), and rewrite the expression as
4737 *
4738 * (-1)^y * (x')^y
4739 *
4740 * We can only calculate (-1)^y for integer values of y, since
4741 * the result is complex if y is not an integer.
4742 */
4743 if (get_neg(ext_))
4744 {
4745 size_t idx;
4746 int units_dig;
4747
4748 /* copy x into r2 */
4749 copy_val(ext2, ext_, FALSE);
4750
4751 /*
4752 * calculate x' = (-x) - since x is negative, this will
4753 * guarantee that x' is positive
4754 */
4755 set_neg(ext2, FALSE);
4756
4757 /*
4758 * make sure y is an integer - start at the first digit
4759 * after the decimal point and check for any non-zero digits
4760 */
4761 idx = (get_exp(val2_ext) < 0 ? 0 : (size_t)get_exp(val2_ext));
4762 for ( ; idx < get_prec(val2_ext) ; ++idx)
4763 {
4764 /* if this digit isn't a zero, it's not an integer */
4765 if (get_dig(val2_ext, idx) != 0)
4766 {
4767 /* y isn't an integer, so we can't calculate (-1)^y */
4768 err_throw(VMERR_OUT_OF_RANGE);
4769 }
4770 }
4771
4772 /* get the first digit to the left of the decimal point */
4773 if (get_exp(val2_ext) <= 0
4774 || (size_t)get_exp(val2_ext) > get_prec(val2_ext))
4775 {
4776 /* the units digit isn't represented - zero is implied */
4777 units_dig = 0;
4778 }
4779 else
4780 {
4781 /* get the digit */
4782 units_dig = get_dig(val2_ext, (size_t)get_exp(val2_ext) - 1);
4783 }
4784
4785 /*
4786 * if the units digit is even, the result will be positive;
4787 * if it's odd, the result will be negative
4788 */
4789 result_neg = ((units_dig & 1) != 0);
4790
4791 /* calculate ln(x') into r1 */
4792 compute_ln_into(vmg_ ext1, ext2);
4793 }
4794 else
4795 {
4796 /* calculate ln(x) into r1 */
4797 compute_ln_into(vmg_ ext1, ext_);
4798
4799 /* the result will be positive */
4800 result_neg = FALSE;
4801 }
4802
4803 /* calculate y * ln(x) into r2 */
4804 compute_prod_into(ext2, val2_ext, ext1);
4805
4806 /* calculate exp(r2) = exp(y*ln(x)) = x^y into r1 */
4807 compute_exp_into(vmg_ ext1, ext2);
4808
4809 /* negate the result if we had a negative x and an odd power */
4810 if (result_neg)
4811 negate(ext1);
4812
4813 /* save the result, rounding as needed */
4814 copy_val(new_ext, ext1, TRUE);
4815 }
4816 err_finally
4817 {
4818 /* release our temporary registers */
4819 release_temp_regs(vmg_ 2, hdl1, hdl2);
4820 }
4821 err_end;
4822
4823 done:
4824 /* discard the GC protection */
4825 G_stk->discard(2);
4826
4827 /* handled */
4828 return TRUE;
4829 }
4830
4831 /*
4832 * property evaluator - sinh
4833 */
getp_sinh(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)4834 int CVmObjBigNum::getp_sinh(VMG_ vm_obj_id_t self,
4835 vm_val_t *retval, uint *argc)
4836 {
4837 /* calculate the hyperbolic sine using the common evaluator */
4838 return calc_sinhcosh(vmg_ self, retval, argc, FALSE, FALSE);
4839 }
4840
4841 /*
4842 * property evaluator - cosh
4843 */
getp_cosh(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)4844 int CVmObjBigNum::getp_cosh(VMG_ vm_obj_id_t self,
4845 vm_val_t *retval, uint *argc)
4846 {
4847 /* calculate the hyperbolic cosine using the common evaluator */
4848 return calc_sinhcosh(vmg_ self, retval, argc, TRUE, FALSE);
4849 }
4850
getp_tanh(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc)4851 int CVmObjBigNum::getp_tanh(VMG_ vm_obj_id_t self,
4852 vm_val_t *retval, uint *argc)
4853 {
4854 /* calculate the hyperbolic tangent using the common evaluator */
4855 return calc_sinhcosh(vmg_ self, retval, argc, FALSE, TRUE);
4856 }
4857
4858 /*
4859 * common property evaluator - compute a hyperbolic sine or cosine, or
4860 * do both to calculate a hyperbolic tangent
4861 */
calc_sinhcosh(VMG_ vm_obj_id_t self,vm_val_t * retval,uint * argc,int is_cosh,int is_tanh)4862 int CVmObjBigNum::calc_sinhcosh(VMG_ vm_obj_id_t self,
4863 vm_val_t *retval, uint *argc,
4864 int is_cosh, int is_tanh)
4865 {
4866 char *new_ext;
4867
4868 /* check arguments and allocate the return value */
4869 if (setup_getp_0(vmg_ self, retval, argc, &new_ext))
4870 return TRUE;
4871
4872 /* calculate the result */
4873 compute_sinhcosh_into(vmg_ new_ext, ext_, is_cosh, is_tanh);
4874
4875 /* discard the GC protection */
4876 G_stk->discard();
4877
4878 /* handled */
4879 return TRUE;
4880 }
4881
4882 /* ------------------------------------------------------------------------ */
4883 /*
4884 * Compute a natural logarithm
4885 */
compute_ln_into(VMG_ char * dst,const char * src)4886 void CVmObjBigNum::compute_ln_into(VMG_ char *dst, const char *src)
4887 {
4888 uint hdl1, hdl2, hdl3, hdl4, hdl5;
4889 char *ext1, *ext2, *ext3, *ext4, *ext5;
4890 size_t prec = get_prec(dst);
4891 const char *ln10;
4892
4893 /* cache the value of ln(10) */
4894 ln10 = cache_ln10(vmg_ prec + 3);
4895
4896 /* if the source value is zero, it's an error */
4897 if (is_zero(src))
4898 err_throw(VMERR_OUT_OF_RANGE);
4899
4900 /* allocate some temporary registers */
4901 alloc_temp_regs(vmg_ prec + 3, 5,
4902 &ext1, &hdl1, &ext2, &hdl2, &ext3, &hdl3,
4903 &ext4, &hdl4, &ext5, &hdl5);
4904
4905 /* catch errors so we can be sure to free registers */
4906 err_try
4907 {
4908 int src_exp;
4909
4910 /*
4911 * Observe that we store our values as x = a*10^b. We can thus
4912 * rewrite ln(x) as ln(a*10^b) = ln(a)+ln(10^b) =
4913 * log(a)+b*ln(10). a is the mantissa, which due to
4914 * normalization is in the range 0.1 <= a < 1.0. Thus, a is an
4915 * ideal argument for the Taylor series. So, we can simply
4916 * compute ln(mantissa), then add it to ln(10)*exponent.
4917 */
4918
4919 /* copy our argument into r1 */
4920 copy_val(ext1, src, FALSE);
4921
4922 /*
4923 * remember the original exponent, and set the exponent of the
4924 * series argument to zero - we'll correct for this later by
4925 * adding the log of the exponent to the series result
4926 */
4927 src_exp = get_exp(ext1);
4928 set_exp(ext1, 0);
4929
4930 /*
4931 * if the lead digit of the mantissa is 1, multiply the number
4932 * by ten and adjust the exponent accordingly - the series is
4933 * especially good at values near 1
4934 */
4935 if (get_dig(ext1, 0) == 1)
4936 {
4937 set_exp(ext1, 1);
4938 --src_exp;
4939 }
4940
4941 /* compute the series expansion */
4942 ext1 = compute_ln_series_into(vmg_ ext1, ext2, ext3, ext4, ext5);
4943
4944 /* add in the input exponent, properly adjusted */
4945 if (src_exp != 0)
4946 {
4947 /* get ln10 into r2 */
4948 copy_val(ext2, ln10, TRUE);
4949
4950 /* apply the exponent's sign if it's negative */
4951 if (src_exp < 0)
4952 {
4953 set_neg(ext2, TRUE);
4954 src_exp = -src_exp;
4955 }
4956
4957 /* multiply by the exponent */
4958 mul_by_long(ext2, src_exp);
4959
4960 /* add this value to the series expansion value */
4961 compute_sum_into(ext3, ext1, ext2);
4962
4963 /* use this as the result */
4964 ext1 = ext3;
4965 }
4966
4967 /* copy the result, rounding if needed */
4968 copy_val(dst, ext1, TRUE);
4969 }
4970 err_finally
4971 {
4972 /* release our temporary registers */
4973 release_temp_regs(vmg_ 5, hdl1, hdl2, hdl3, hdl4, hdl5);
4974 }
4975 err_end;
4976 }
4977
4978 /*
4979 * Compute the natural log series. The argument value, initially in
4980 * ext1, should be adjusted to a small value before this is called to
4981 * ensure quick series convergence.
4982 */
compute_ln_series_into(VMG_ char * ext1,char * ext2,char * ext3,char * ext4,char * ext5)4983 char *CVmObjBigNum::compute_ln_series_into(VMG_ char *ext1, char *ext2,
4984 char *ext3, char *ext4,
4985 char *ext5)
4986 {
4987 ulong n;
4988 char *tmp;
4989
4990 /* start at the first term of the series */
4991 n = 1;
4992
4993 /* subtract one from r1 to yield (x-1) in r2 */
4994 compute_abs_diff_into(ext2, ext1, get_one());
4995
4996 /* add one to r1 to yield (x+1) */
4997 increment_abs(ext1);
4998
4999 /* compute (x-1)/(x+1) into r3 - this will be our current power */
5000 compute_quotient_into(vmg_ ext3, 0, ext2, ext1);
5001
5002 /*
5003 * compute ((x-1)/(x+1))^2 into r4 - we'll multiply r3 by this on
5004 * each iteration to produce the next required power, since we only
5005 * need the odd powers
5006 */
5007 compute_prod_into(ext4, ext3, ext3);
5008
5009 /* start out with the first power in our accumulator (r1) */
5010 copy_val(ext1, ext3, FALSE);
5011
5012 /* iterate until we have a good enough answer */
5013 for (;;)
5014 {
5015 /* compute the next power into r2 */
5016 compute_prod_into(ext2, ext3, ext4);
5017
5018 /* copy the result into our current power register, r3 */
5019 copy_val(ext3, ext2, FALSE);
5020
5021 /* advance n */
5022 n += 2;
5023
5024 /* divide the power by n */
5025 div_by_long(ext2, n);
5026
5027 /* if it's too small to notice, we're done */
5028 if (is_zero(ext2)
5029 || get_exp(ext1) - get_exp(ext2) > (int)get_prec(ext1))
5030 break;
5031
5032 /* compute the sum with our accumulator into r5 */
5033 compute_sum_into(ext5, ext1, ext2);
5034
5035 /* swap r5 and r1 - the new sum is our new accumulator */
5036 tmp = ext5;
5037 ext5 = ext1;
5038 ext1 = tmp;
5039 }
5040
5041 /* multiply the result of the series by 2 */
5042 mul_by_long(ext1, 2);
5043
5044 /* return the register containing the result */
5045 return ext1;
5046 }
5047
5048 /*
5049 * Compute e^x, where e is the base of the natural logarithm (2.818...)
5050 */
compute_exp_into(VMG_ char * dst,const char * src)5051 void CVmObjBigNum::compute_exp_into(VMG_ char *dst, const char *src)
5052 {
5053 uint hdl1, hdl2, hdl3, hdl4, hdl5, hdl6;
5054 char *ext1, *ext2, *ext3, *ext4, *ext5, *ext6;
5055 size_t prec = get_prec(dst);
5056 const char *ln10;
5057
5058 /* get the constant value of ln10 to the required precision */
5059 ln10 = cache_ln10(vmg_ prec + 3);
5060
5061 /* allocate some temporary registers */
5062 alloc_temp_regs(vmg_ prec + 3, 6,
5063 &ext1, &hdl1, &ext2, &hdl2, &ext3, &hdl3,
5064 &ext4, &hdl4, &ext5, &hdl5, &ext6, &hdl6);
5065
5066 /* catch errors so we can be sure to free registers */
5067 err_try
5068 {
5069 char *tmp;
5070 ulong n;
5071 ulong n_fact;
5072 int use_int_fact;
5073 long new_exp;
5074 size_t idx;
5075 size_t decpt_idx;
5076 int twos;
5077
5078 /*
5079 * We want to calculate e^x. Observe that a^x = e^(x*ln(x)), so
5080 * 10^x = e^(x*ln(10)). We can rewrite our desired expression
5081 * e^x as
5082 *
5083 * e^[(x/ln(10)) * ln(10)]
5084 *
5085 * which is thus
5086 *
5087 * 10^(x / ln(10))
5088 *
5089 * We store our numbers as a mantissa times 10 raised to an
5090 * integer exponent. Clearly the exponent of 10 in the formula
5091 * above is not always an integer (in fact, it's extremely rare
5092 * that it would be an integer), so we still have more work to
5093 * do. What we must do is obtain an integer exponent. So, let
5094 * us define y = x/ln(10); we can now rewrite the above as
5095 *
5096 * 10^(int(y) + frac(y))
5097 *
5098 * where int(y) is the integer portion of y and frac(y) is the
5099 * fractional portion (y - int(y)). We can rewrite the above as
5100 *
5101 * 10^frac(y) * 10^int(y)
5102 *
5103 * which we can further rewrite as
5104 *
5105 * e^(frac(y)*ln(10)) * 10^int(y)
5106 *
5107 * We haven't made the problem of finding an exponential
5108 * disappear, but we've reduced the argument to a very
5109 * manageable range, which is important because it makes the
5110 * Taylor series converge quickly. Furthermore, it's extremely
5111 * inexpensive to separate out the problem like this, since it
5112 * falls quite naturally out of the representation we use, so it
5113 * doesn't add much overhead to do this preparation work.
5114 */
5115
5116 /* first, calculate x/ln(10) into r1 */
5117 compute_quotient_into(vmg_ ext1, 0, src, ln10);
5118
5119 /*
5120 * compute the integer portion of x/ln(10) - it has to fit in a
5121 * 16-bit integer, (-32768 to +32767), because this is going to
5122 * be the exponent of the result (or roughly so, anyway)
5123 */
5124 decpt_idx = get_exp(ext1) >= 0 ? (size_t)get_exp(ext1) : 0;
5125 for (new_exp = 0, idx = 0 ; idx < decpt_idx ; ++idx)
5126 {
5127 int dig;
5128
5129 /*
5130 * get this digit if it's represented; if not, it's an
5131 * implied trailing zero
5132 */
5133 dig = (idx < get_prec(ext1) ? get_dig(ext1, idx) : 0);
5134
5135 /* add this digit into the accumulator */
5136 new_exp *= 10;
5137 new_exp += dig;
5138
5139 /*
5140 * Make sure we're still in range. Note that, because our
5141 * representation is 0.dddd*10^x, we need one more factor of
5142 * ten than you'd think here, the adjust of the range from
5143 * the expected -32768..32767
5144 */
5145 if (new_exp > (get_neg(ext1) ? 32769L : 32766L))
5146 err_throw(VMERR_NUM_OVERFLOW);
5147
5148 /*
5149 * zero out this digit, so that when we're done r1 has the
5150 * fractional part only
5151 */
5152 if (idx < get_prec(ext1))
5153 set_dig(ext1, idx, 0);
5154 }
5155
5156 /* negate the exponent value if the source value is negative */
5157 if (get_neg(ext1))
5158 new_exp = -new_exp;
5159
5160 /* normalize the fractional part, which remains in ext1 */
5161 normalize(ext1);
5162
5163 /*
5164 * Multiply it by ln10, storing the result in r3. This is the
5165 * value we will use with the Taylor series.
5166 */
5167 compute_prod_into(ext3, ext1, ln10);
5168
5169 /*
5170 * While our input value is greater than 0.5, divide it by two
5171 * to make it smaller than 0.5. This will speed up the series
5172 * convergence. When we're done, we'll correct for the
5173 * divisions my squaring the result: e^2x = (e^x)^2
5174 */
5175 copy_val(ext1, get_one(), FALSE);
5176 div_by_long(ext1, 2);
5177 for (twos = 0 ; compare_abs(ext3, ext1) > 0 ; ++twos)
5178 div_by_long(ext3, 2);
5179
5180 /* start with 1 in our accumulator (r1) */
5181 copy_val(ext1, get_one(), FALSE);
5182
5183 /* copy our series argument into the current-power register (r2) */
5184 copy_val(ext2, ext3, FALSE);
5185
5186 /* start with 1 in our factorial register (r4) */
5187 copy_val(ext4, get_one(), FALSE);
5188
5189 /* start at the first term */
5190 n = 1;
5191 n_fact = 1;
5192 use_int_fact = TRUE;
5193
5194 /* go until we reach the required precision */
5195 for (;;)
5196 {
5197 /* for efficiency, try integer division */
5198 if (use_int_fact)
5199 {
5200 /*
5201 * we can still fit the factorial in an integer - divide
5202 * by the integer value of n!, since it's a lot faster
5203 */
5204 copy_val(ext5, ext2, FALSE);
5205 div_by_long(ext5, n_fact);
5206
5207 /* calculate the next n! integer, if it'll fit in a long */
5208 if (n_fact > LONG_MAX/(n+1))
5209 {
5210 /*
5211 * it'll be too big next time - we'll have to start
5212 * using the full quotient calculation
5213 */
5214 use_int_fact = FALSE;
5215 }
5216 else
5217 {
5218 /* it'll still fit - calculate the next n! */
5219 n_fact *= (n+1);
5220 }
5221 }
5222 else
5223 {
5224 /* compute x^n/n! (r2/r4) into r5 */
5225 compute_quotient_into(vmg_ ext5, 0, ext2, ext4);
5226 }
5227
5228 /* if we're below the required precision, we're done */
5229 if (is_zero(ext5)
5230 || get_exp(ext1) - get_exp(ext5) > (int)get_prec(ext1))
5231 break;
5232
5233 /* compute the sum of the accumulator and this term into r6 */
5234 compute_sum_into(ext6, ext1, ext5);
5235
5236 /* swap the result into the accumulator */
5237 tmp = ext1;
5238 ext1 = ext6;
5239 ext6 = tmp;
5240
5241 /* on to the next term */
5242 ++n;
5243
5244 /* compute the next factorial value */
5245 mul_by_long(ext4, n);
5246
5247 /* compute the next power of x' into r5 */
5248 compute_prod_into(ext5, ext2, ext3);
5249
5250 /* swap the result into our current-power register (r2) */
5251 tmp = ext2;
5252 ext2 = ext5;
5253 ext5 = tmp;
5254 }
5255
5256 /* square the result as many times as we halved the argument */
5257 for ( ; twos != 0 ; --twos)
5258 {
5259 /* compute the square of r1 into r2 */
5260 compute_prod_into(ext2, ext1, ext1);
5261
5262 /* swap the result back into r1 */
5263 tmp = ext1;
5264 ext1 = ext2;
5265 ext2 = tmp;
5266 }
5267
5268 /*
5269 * set up our 10's exponent value - this is simply 1*10^new_exp,
5270 * which we calculated earlier (which we represent as
5271 * 0.1*10^(new_exp+1)
5272 */
5273 copy_val(ext2, get_one(), FALSE);
5274 set_exp(ext2, (int)(new_exp + 1));
5275
5276 /* multiply by the 10's exponent value */
5277 compute_prod_into(ext3, ext1, ext2);
5278
5279 /* copy the result into the output register, rounding as needed */
5280 copy_val(dst, ext3, TRUE);
5281 }
5282 err_finally
5283 {
5284 /* release our temporary registers */
5285 release_temp_regs(vmg_ 6, hdl1, hdl2, hdl3, hdl4, hdl5, hdl6);
5286 }
5287 err_end;
5288 }
5289
5290 /* ------------------------------------------------------------------------ */
5291 /*
5292 * Compute a hyperbolic sine or cosine
5293 */
compute_sinhcosh_into(VMG_ char * dst,const char * src,int is_cosh,int is_tanh)5294 void CVmObjBigNum::compute_sinhcosh_into(VMG_ char *dst, const char *src,
5295 int is_cosh, int is_tanh)
5296 {
5297 size_t prec = get_prec(dst);
5298 uint hdl1, hdl2, hdl3, hdl4;
5299 char *ext1, *ext2, *ext3, *ext4;
5300
5301 /* allocate some temporary registers */
5302 alloc_temp_regs(vmg_ prec + 3, 4,
5303 &ext1, &hdl1, &ext2, &hdl2, &ext3, &hdl3,
5304 &ext4, &hdl4);
5305
5306 /* catch errors so we can be sure to free registers */
5307 err_try
5308 {
5309 /*
5310 * sinh(x) = (e^x - e^(-x))/2
5311 *
5312 * cosh(x) = (e^x + e^(-x))/2
5313 *
5314 * The two differ only in the sign of the e^(-x) term, so most
5315 * of the calculation is the same for both.
5316 */
5317
5318 /* compute e^x into r1 */
5319 compute_exp_into(vmg_ ext1, src);
5320
5321 /*
5322 * rather than calculating e^-x separately, simply invert our
5323 * e^x value to yield e^-x (a simple division is much quicker
5324 * than calculating another exponent, which involves an entire
5325 * taylor series expansion)
5326 */
5327 copy_val(ext3, get_one(), FALSE);
5328 compute_quotient_into(vmg_ ext2, 0, ext3, ext1);
5329
5330 /*
5331 * if we're calculating the tanh, we'll want both the sinh and
5332 * cosh values
5333 */
5334 if (is_tanh)
5335 {
5336 /* add the terms to get the cosh */
5337 compute_sum_into(ext4, ext1, ext2);
5338
5339 /* subtract ext2 to get the sinh */
5340 negate(ext2);
5341 compute_sum_into(ext3, ext1, ext2);
5342
5343 /* tanh is the quotient of sinh/cosh */
5344 compute_quotient_into(vmg_ ext1, 0, ext3, ext4);
5345
5346 /* our result is in ext1 - set ext3 to point there */
5347 ext3 = ext1;
5348 }
5349 else
5350 {
5351 /*
5352 * if this is sinh, the e^-x term is subtracted; if it's
5353 * cosh, it's added
5354 */
5355 if (!is_cosh)
5356 negate(ext2);
5357
5358 /* compute r1 + r2 into r3 (e^x +/- e^(-x)) */
5359 compute_sum_into(ext3, ext1, ext2);
5360
5361 /* halve the result */
5362 div_by_long(ext3, 2);
5363 }
5364
5365 /* store the result */
5366 copy_val(dst, ext3, TRUE);
5367 }
5368 err_finally
5369 {
5370 /* release our temporary registers */
5371 release_temp_regs(vmg_ 4, hdl1, hdl2, hdl3, hdl4);
5372 }
5373 err_end;
5374 }
5375
5376 /* ------------------------------------------------------------------------ */
5377 /*
5378 * Cache the natural logarithm of 10 to the given precision and return
5379 * the value
5380 */
cache_ln10(VMG_ size_t prec)5381 const char *CVmObjBigNum::cache_ln10(VMG_ size_t prec)
5382 {
5383 char *ext;
5384 int expanded;
5385 uint hdl1, hdl2, hdl3, hdl4, hdl5;
5386 char *ext1, *ext2, *ext3, *ext4, *ext5;
5387 static const unsigned char ten[] =
5388 {
5389 /* number of digits */
5390 0x01, 0x00,
5391
5392 /* exponent (0.1 * 10^2 = 10) */
5393 0x02, 0x00,
5394
5395 /* flags */
5396 0x00,
5397
5398 /* mantissa */
5399 0x10
5400 };
5401
5402 /*
5403 * round up the precision a bit to ensure that we don't have to
5404 * repeatedly recalculate this value if we're asked for a cluster of
5405 * similar precisions
5406 */
5407 prec = (prec + 7) & ~7;
5408
5409 /* get the ln10 cache register */
5410 ext = G_bignum_cache->get_ln10_reg(calc_alloc(prec), &expanded);
5411
5412 /* if that failed, throw an error */
5413 if (ext == 0)
5414 err_throw(VMERR_OUT_OF_MEMORY);
5415
5416 /*
5417 * if we didn't have to allocate more memory, and the value in the
5418 * register has at least the required precision, return the cached
5419 * value
5420 */
5421 if (!expanded && get_prec(ext) >= prec)
5422 return ext;
5423
5424 /*
5425 * we have reallocated the register, or we just didn't have enough
5426 * precision in the old value - set the new precision
5427 */
5428 set_prec(ext, prec);
5429
5430 /* allocate some temporary registers */
5431 alloc_temp_regs(vmg_ prec + 3, 5,
5432 &ext1, &hdl1, &ext2, &hdl2, &ext3, &hdl3,
5433 &ext4, &hdl4, &ext5, &hdl5);
5434
5435 /* catch errors so we can be sure to free registers */
5436 err_try
5437 {
5438 /*
5439 * Compute sqrt(10) - 10 is too large for the series to
5440 * converge, but sqrt(10) is good. We'll correct for this later
5441 * by doubling the result of the series expansion, which gives
5442 * us the correct result: ln(a^b) = b*ln(a), and sqrt(x) =
5443 * x^(1/2), hence ln(sqrt(x)) = ln(x)/2, which means that ln(x)
5444 * = 2*ln(sqrt(x)).
5445 */
5446
5447 /* compute sqrt(10), for quick series convergence */
5448 copy_val(ext2, (const char *)ten, FALSE);
5449 compute_sqrt_into(vmg_ ext1, ext2);
5450
5451 /* compute the series expansion */
5452 ext1 = compute_ln_series_into(vmg_ ext1, ext2, ext3, ext4, ext5);
5453
5454 /* double the result (to adjust for the sqrt) */
5455 mul_by_long(ext1, 2);
5456
5457 /* store the result in the cache */
5458 copy_val(ext, ext1, TRUE);
5459 }
5460 err_finally
5461 {
5462 /* release our temporary registers */
5463 release_temp_regs(vmg_ 5, hdl1, hdl2, hdl3, hdl4, hdl5);
5464 }
5465 err_end;
5466
5467 /* return the register pointer */
5468 return ext;
5469 }
5470
5471 /* ------------------------------------------------------------------------ */
5472 /*
5473 * Cache pi to the required precision
5474 */
cache_pi(VMG_ size_t prec)5475 const char *CVmObjBigNum::cache_pi(VMG_ size_t prec)
5476 {
5477 char *ext;
5478 int expanded;
5479 uint hdl1, hdl2, hdl3, hdl4, hdl5;
5480 char *ext1, *ext2, *ext3, *ext4, *ext5;
5481
5482 /*
5483 * round up the precision a bit to ensure that we don't have to
5484 * repeatedly recalculate this value if we're asked for a cluster of
5485 * similar precisions
5486 */
5487 prec = (prec + 7) & ~7;
5488
5489 /* get the pi cache register */
5490 ext = G_bignum_cache->get_pi_reg(calc_alloc(prec), &expanded);
5491
5492 /* if that failed, throw an error */
5493 if (ext == 0)
5494 err_throw(VMERR_OUT_OF_MEMORY);
5495
5496 /*
5497 * if we didn't have to allocate more memory, and the value in the
5498 * register has at least the required precision, return the cached
5499 * value
5500 */
5501 if (!expanded && get_prec(ext) >= prec)
5502 return ext;
5503
5504 /*
5505 * we have reallocated the register, or we just didn't have enough
5506 * precision in the old value - set the new precision
5507 */
5508 set_prec(ext, prec);
5509
5510 /* allocate some temporary registers */
5511 alloc_temp_regs(vmg_ prec + 3, 5,
5512 &ext1, &hdl1, &ext2, &hdl2, &ext3, &hdl3,
5513 &ext4, &hdl4, &ext5, &hdl5);
5514
5515 /* catch errors so we can be sure to free registers */
5516 err_try
5517 {
5518 /*
5519 * Compute the arctangent of 1. To do this, rewrite arctan(x) =
5520 * arccos(1/sqrt(1+x^2)), so x=1 gives us arccos(1/sqrt(2)).
5521 *
5522 * Now, arccos(x) = pi/2 - arcsin(x), but this would defeat the
5523 * purpose since we want to calculate pi, thus we don't want to
5524 * depend upon pi in our calculations. Fortunately, arcsin(x) =
5525 * pi/2 - arcsin(sqrt(1-x^2)), hence
5526 *
5527 * pi/2 - arcsin(x) = pi/2 - (pi/2 - arcsin(sqrt(1-x^2)))
5528 *. = arcsin(sqrt(1-x^2))
5529 *
5530 * So, we can rewrite arccos(1/sqrt(2)) as arcsin(sqrt(1/2)).
5531 */
5532
5533 /* compute 1/2 into r2 */
5534 copy_val(ext2, get_one(), FALSE);
5535 div_by_long(ext2, 2);
5536
5537 /* compute sqrt(1/2) into r1 */
5538 compute_sqrt_into(vmg_ ext1, ext2);
5539
5540 /* calculate the arcsin series for sqrt(1/2) */
5541 ext1 = calc_asin_series(ext1, ext2, ext3, ext4, ext5);
5542
5543 /* multiply the result by 4 */
5544 mul_by_long(ext1, 4);
5545
5546 /* store the result in the cache */
5547 copy_val(ext, ext1, TRUE);
5548 }
5549 err_finally
5550 {
5551 /* release our temporary registers */
5552 release_temp_regs(vmg_ 5, hdl1, hdl2, hdl3, hdl4, hdl5);
5553 }
5554 err_end;
5555
5556 /* return the register pointer */
5557 return ext;
5558 }
5559
5560 /* ------------------------------------------------------------------------ */
5561 /*
5562 * Cache e to the required precision
5563 */
cache_e(VMG_ size_t prec)5564 const char *CVmObjBigNum::cache_e(VMG_ size_t prec)
5565 {
5566 char *ext;
5567 int expanded;
5568
5569 /*
5570 * round up the precision a bit to ensure that we don't have to
5571 * repeatedly recalculate this value if we're asked for a cluster of
5572 * similar precisions
5573 */
5574 prec = (prec + 7) & ~7;
5575
5576 /* get the e cache register */
5577 ext = G_bignum_cache->get_e_reg(calc_alloc(prec), &expanded);
5578
5579 /* if that failed, throw an error */
5580 if (ext == 0)
5581 err_throw(VMERR_OUT_OF_MEMORY);
5582
5583 /*
5584 * if we didn't have to allocate more memory, and the value in the
5585 * register has at least the required precision, return the cached
5586 * value
5587 */
5588 if (!expanded && get_prec(ext) >= prec)
5589 return ext;
5590
5591 /*
5592 * we have reallocated the register, or we just didn't have enough
5593 * precision in the old value - set the new precision
5594 */
5595 set_prec(ext, prec);
5596
5597 /* exponentiate the base of the natural logarithm to the power 1 */
5598 compute_exp_into(vmg_ ext, get_one());
5599
5600 /* return the register pointer */
5601 return ext;
5602 }
5603
5604 /* ------------------------------------------------------------------------ */
5605 /*
5606 * Compute a square root. We use Newton's method (a reliable old method
5607 * for extracting square roots, made better by the fact that, unlike the
5608 * method's inventor, we are fortunate to have an electronic computer to
5609 * carry out the tedious parts of the calculation for us).
5610 */
compute_sqrt_into(VMG_ char * dst,const char * src)5611 void CVmObjBigNum::compute_sqrt_into(VMG_ char *dst, const char *src)
5612 {
5613 uint hdl1, hdl2, hdl3, hdl4;
5614 char *ext1, *ext2, *ext3, *ext4;
5615 size_t dst_prec = get_prec(dst);
5616
5617 /* if the value is negative, it's an error */
5618 if (get_neg(src))
5619 err_throw(VMERR_OUT_OF_RANGE);
5620
5621 /* allocate our scratchpad registers */
5622 alloc_temp_regs(vmg_ get_prec(dst) + 3, 4,
5623 &ext1, &hdl1, &ext2, &hdl2, &ext3, &hdl3,
5624 &ext4, &hdl4);
5625
5626 /* catch errors so we can free our registers */
5627 err_try
5628 {
5629 /*
5630 * Compute our initial guess. Since our number is represented as
5631 *
5632 * (n * 10^exp),
5633 *
5634 * the square root can be written as
5635 *
5636 * sqrt(n) * sqrt(10^exp) = sqrt(n) * 10^(exp/2)
5637 *
5638 * Approximate sqrt(n) as simply n for our initial guess. This
5639 * will get us to the right order of magnitude plus or minus
5640 * one, so we should converge pretty quickly.
5641 *
5642 * If we have an odd exponent, round up and divide the mantissa
5643 * by 2 - this will be something like 0.456e7, which can be
5644 * written as 4.56e6, whose square root is about 2e3, or .2e4.
5645 *
5646 * If we have an even exponent, multiply the mantissa by 2.
5647 * This will be something like .456e8, whose square root is
5648 * about .67e4.
5649 *
5650 * Note that it's well worth the trouble to make a good initial
5651 * approximation, even with the multiply/divide, because these
5652 * operations with longs are much more efficient than the full
5653 * divide we will have to do on each iteration.
5654 */
5655 copy_val(ext1, src, TRUE);
5656 if ((get_exp(ext1) & 1) != 0)
5657 {
5658 /* odd exponent - round up and divide the mantissa by two */
5659 set_exp(ext1, (get_exp(ext1) + 1)/2);
5660 div_by_long(ext1, 2);
5661 }
5662 else
5663 {
5664 /* even exponent - multiply mantissa by two */
5665 set_exp(ext1, get_exp(ext1)/2);
5666 mul_by_long(ext1, 2);
5667 }
5668
5669 /* iterate until we get close enough to the solution */
5670 for (;;)
5671 {
5672 char *tmp;
5673 size_t idx;
5674
5675 /*
5676 * Calculate the next iteration's approximation, noting that r1
5677 * contains the current iteration's value p:
5678 *
5679 * p' = p/2 + src/2p = (p + src/p)/2
5680 *
5681 * Note that if p == 0, we can't compute src/p, so we can't
5682 * iterate any further.
5683 */
5684 if (is_zero(ext1))
5685 break;
5686
5687 /* calculate src/p into r3 */
5688 compute_quotient_into(vmg_ ext3, 0, src, ext1);
5689
5690 /* compute p + src/p into r4 */
5691 compute_sum_into(ext4, ext1, ext3);
5692
5693 /* compute (p + src/p)/2 into r4 */
5694 div_by_long(ext4, 2);
5695
5696 /*
5697 * check for convergence - if the new value equals the old
5698 * value to the precision requested for the result, we are
5699 * at the limit of our ability to distinguish differences in
5700 * future terms, so we can stop
5701 */
5702 if (get_neg(ext1) == get_neg(ext4)
5703 && get_exp(ext1) == get_exp(ext4))
5704 {
5705 /*
5706 * they're the same sign and magnitude - compare the
5707 * digits to see where they first differ
5708 */
5709 for (idx = 0 ; idx < dst_prec + 1 ; ++idx)
5710 {
5711 /* if they differ here, stop scanning */
5712 if (get_dig(ext1, idx) != get_dig(ext4, idx))
5713 break;
5714 }
5715
5716 /*
5717 * if we didn't find any difference up to the output
5718 * precision plus one digit (for rounding), further
5719 * iteration will be of no value
5720 */
5721 if (idx >= dst_prec + 1)
5722 break;
5723 }
5724
5725 /* swap the new value into r1 for the next round */
5726 tmp = ext1;
5727 ext1 = ext4;
5728 ext4 = tmp;
5729 }
5730
5731 /*
5732 * copy the last iteration's value into the destination,
5733 * rounding as needed
5734 */
5735 copy_val(dst, ext1, TRUE);
5736 }
5737 err_finally
5738 {
5739 /* release our temporary registers */
5740 release_temp_regs(vmg_ 4, hdl1, hdl2, hdl3, hdl4);
5741 }
5742 err_end;
5743 }
5744
5745 /* ------------------------------------------------------------------------ */
5746 /*
5747 * Calculate a Taylor expansion for sin(). The argument is in ext1, and
5748 * we'll store the result in new_ext. ext1 through ext7 are temporary
5749 * registers that we'll overwrite with intermediate values.
5750 *
5751 * Before calling this function, the caller should reduce the value in
5752 * ext1 to the range -pi/4 <= ext1 <= pi/4 to ensure quick convergence
5753 * of the series.
5754 */
calc_sin_series(VMG_ char * new_ext,char * ext1,char * ext2,char * ext3,char * ext4,char * ext5,char * ext6,char * ext7)5755 void CVmObjBigNum::calc_sin_series(VMG_ char *new_ext, char *ext1,
5756 char *ext2, char *ext3, char *ext4,
5757 char *ext5, char *ext6, char *ext7)
5758 {
5759 ulong n;
5760 int neg_term;
5761 char *tmp;
5762
5763 /* start with 1! (i.e., 1) in r3 */
5764 copy_val(ext3, get_one(), FALSE);
5765
5766 /*
5767 * calculate x^2 into r7 - we need x, x^3, x^5, etc, so we can just
5768 * multiply the last value by x^2 to get the next power we need
5769 */
5770 compute_prod_into(ext7, ext1, ext1);
5771
5772 /* start with x (reduced mod 2pi) in r5 (our accumulator) */
5773 copy_val(ext5, ext1, FALSE);
5774
5775 /* start at term n=1 in our expansion */
5776 n = 1;
5777
5778 /* the first term is positive */
5779 neg_term = FALSE;
5780
5781 /* go until we have a precise enough value */
5782 for (;;)
5783 {
5784 /*
5785 * move on to the next term: multiply r1 by x^2 into r2 to yield
5786 * the next power of x that we require
5787 */
5788 compute_prod_into(ext2, ext1, ext7);
5789
5790 /* swap r1 and r2 - r2 has the next power we need now */
5791 tmp = ext1;
5792 ext1 = ext2;
5793 ext2 = tmp;
5794
5795 /*
5796 * multiply r3 by (n+1)*(n+2) - it currently has n!, so this
5797 * will yield (n+2)!
5798 */
5799 mul_by_long(ext3, (n+1)*(n+2));
5800
5801 /* advance n to the next term */
5802 n += 2;
5803
5804 /* each term swaps signs */
5805 neg_term = !neg_term;
5806
5807 /* divide r1 by r3 to yield the next term in our series in r4 */
5808 compute_quotient_into(vmg_ ext4, 0, ext1, ext3);
5809
5810 /*
5811 * if this value is too small to count in our accumulator, we're
5812 * done
5813 */
5814 if (is_zero(ext4)
5815 || get_exp(ext5) - get_exp(ext4) > (int)get_prec(ext4))
5816 break;
5817
5818 /*
5819 * invert the sign of the term in r4 if this is a negative term
5820 */
5821 if (neg_term)
5822 set_neg(ext4, !get_neg(ext4));
5823
5824 /* add r4 to our running accumulator in r5, yielding r6 */
5825 compute_sum_into(ext6, ext5, ext4);
5826
5827 /* swap r5 and r6 to put the accumulated value back into r5 */
5828 tmp = ext5;
5829 ext5 = ext6;
5830 ext6 = tmp;
5831 }
5832
5833 /* we're done - store the result, rounding if necessary */
5834 copy_val(new_ext, ext5, TRUE);
5835 }
5836
5837 /*
5838 * Calculate a Taylor expansion for cos(). The argument is in ext1, and
5839 * we'll store the result in new_ext. ext1 through ext7 are temporary
5840 * registers that we'll overwrite with intermediate values.
5841 *
5842 * Before calling this function, the caller should reduce the value in
5843 * ext1 to the range -pi/4 <= ext1 <= pi/4 to ensure quick convergence
5844 * of the series.
5845 */
calc_cos_series(VMG_ char * new_ext,char * ext1,char * ext2,char * ext3,char * ext4,char * ext5,char * ext6,char * ext7)5846 void CVmObjBigNum::calc_cos_series(VMG_ char *new_ext, char *ext1,
5847 char *ext2, char *ext3, char *ext4,
5848 char *ext5, char *ext6, char *ext7)
5849 {
5850 ulong n;
5851 int neg_term;
5852 char *tmp;
5853
5854 /* start with 2! (i.e., 2) in r3 */
5855 copy_val(ext3, get_one(), FALSE);
5856 set_dig(ext3, 0, 2);
5857
5858 /*
5859 * calculate x^2 into r7 - we need x^2, x^4, x^6, etc, so we can
5860 * just multiply the last value by x^2 to get the next power we need
5861 */
5862 compute_prod_into(ext7, ext1, ext1);
5863
5864 /*
5865 * the first power we need is x^2, so copy it into r1 (our current
5866 * power of x register)
5867 */
5868 copy_val(ext1, ext7, FALSE);
5869
5870 /*
5871 * start with 1 in r5, the accumulator (the first term of the series
5872 * is the constant value 1)
5873 */
5874 copy_val(ext5, get_one(), FALSE);
5875
5876 /*
5877 * start at term n=2 in our expansion (the first term is just the
5878 * constant value 1, so we can start at the second term)
5879 */
5880 n = 2;
5881
5882 /*
5883 * the first term we calculate (i.e., the second term in the actual
5884 * series) is negative
5885 */
5886 neg_term = TRUE;
5887
5888 /* go until we have a precise enough value */
5889 for (;;)
5890 {
5891 /* divide r1 by r3 to yield the next term in our series in r4 */
5892 compute_quotient_into(vmg_ ext4, 0, ext1, ext3);
5893
5894 /*
5895 * if this value is too small to count in our accumulator, we're
5896 * done
5897 */
5898 if (is_zero(ext4)
5899 || get_exp(ext5) - get_exp(ext4) > (int)get_prec(ext4))
5900 break;
5901
5902 /*
5903 * invert the sign of the term in r4 if this is a negative term
5904 */
5905 if (neg_term)
5906 set_neg(ext4, !get_neg(ext4));
5907
5908 /* add r4 to our running accumulator in r5, yielding r6 */
5909 compute_sum_into(ext6, ext5, ext4);
5910
5911 /* swap r5 and r6 to put the accumulated value back into r5 */
5912 tmp = ext5;
5913 ext5 = ext6;
5914 ext6 = tmp;
5915
5916 /*
5917 * move on to the next term: multiply r1 by x^2 into r2 to yield
5918 * the next power of x that we require
5919 */
5920 compute_prod_into(ext2, ext1, ext7);
5921
5922 /* swap r1 and r2 to put our next required power back in r1 */
5923 tmp = ext1;
5924 ext1 = ext2;
5925 ext2 = tmp;
5926
5927 /*
5928 * multiply r3 by (n+1)*(n+2) - it currently has n!, so this
5929 * will yield (n+2)!
5930 */
5931 mul_by_long(ext3, (n+1)*(n+2));
5932
5933 /* advance n to the next term */
5934 n += 2;
5935
5936 /* each term swaps signs */
5937 neg_term = !neg_term;
5938 }
5939
5940 /* we're done - store the result, rounding if necessary */
5941 copy_val(new_ext, ext5, TRUE);
5942 }
5943
5944 /* ------------------------------------------------------------------------ */
5945 /*
5946 * add a value
5947 */
add_val(VMG_ vm_val_t * result,vm_obj_id_t self,const vm_val_t * val)5948 void CVmObjBigNum::add_val(VMG_ vm_val_t *result,
5949 vm_obj_id_t self, const vm_val_t *val)
5950 {
5951 vm_val_t val2;
5952
5953 /* convert it */
5954 val2 = *val;
5955 if (!cvt_to_bignum(vmg_ self, &val2))
5956 {
5957 /* this type is not convertible to BigNumber, so we can't add it */
5958 err_throw(VMERR_BAD_TYPE_ADD);
5959 }
5960
5961 /* push 'self' and the other value to protect against GC */
5962 G_stk->push()->set_obj(self);
5963 G_stk->push(&val2);
5964
5965 /* compute the sum */
5966 compute_sum(vmg_ result, ext_, get_objid_ext(vmg_ val2.val.obj));
5967
5968 /* discard the GC protection items */
5969 G_stk->discard(2);
5970 }
5971
5972 /*
5973 * subtract a value
5974 */
sub_val(VMG_ vm_val_t * result,vm_obj_id_t self,const vm_val_t * val)5975 void CVmObjBigNum::sub_val(VMG_ vm_val_t *result,
5976 vm_obj_id_t self, const vm_val_t *val)
5977 {
5978 vm_val_t val2;
5979
5980 /* convert it */
5981 val2 = *val;
5982 if (!cvt_to_bignum(vmg_ self, &val2))
5983 {
5984 /* this type is not convertible to BigNumber, so we can't use it */
5985 err_throw(VMERR_BAD_TYPE_SUB);
5986 }
5987
5988 /* push 'self' and the other value to protect against GC */
5989 G_stk->push()->set_obj(self);
5990 G_stk->push(&val2);
5991
5992 /* compute the difference */
5993 compute_diff(vmg_ result, ext_, get_objid_ext(vmg_ val2.val.obj));
5994
5995 /* discard the GC protection items */
5996 G_stk->discard(2);
5997 }
5998
5999 /*
6000 * multiply a value
6001 */
mul_val(VMG_ vm_val_t * result,vm_obj_id_t self,const vm_val_t * val)6002 void CVmObjBigNum::mul_val(VMG_ vm_val_t *result,
6003 vm_obj_id_t self, const vm_val_t *val)
6004 {
6005 vm_val_t val2;
6006
6007 /* convert it */
6008 val2 = *val;
6009 if (!cvt_to_bignum(vmg_ self, &val2))
6010 {
6011 /* this type is not convertible to BigNumber, so we can't add it */
6012 err_throw(VMERR_BAD_TYPE_ADD);
6013 }
6014
6015 /* push 'self' and the other value to protect against GC */
6016 G_stk->push()->set_obj(self);
6017 G_stk->push(&val2);
6018
6019 /* compute the product */
6020 compute_prod(vmg_ result, ext_, get_objid_ext(vmg_ val2.val.obj));
6021
6022 /* discard the GC protection items */
6023 G_stk->discard(2);
6024 }
6025
6026 /*
6027 * divide a value
6028 */
div_val(VMG_ vm_val_t * result,vm_obj_id_t self,const vm_val_t * val)6029 void CVmObjBigNum::div_val(VMG_ vm_val_t *result,
6030 vm_obj_id_t self, const vm_val_t *val)
6031 {
6032 vm_val_t val2;
6033
6034 /* convert it */
6035 val2 = *val;
6036 if (!cvt_to_bignum(vmg_ self, &val2))
6037 {
6038 /* this type is not convertible to BigNumber, so we can't add it */
6039 err_throw(VMERR_BAD_TYPE_ADD);
6040 }
6041
6042 /* push 'self' and the other value to protect against GC */
6043 G_stk->push()->set_obj(self);
6044 G_stk->push(&val2);
6045
6046 /* compute the quotient */
6047 compute_quotient(vmg_ result, ext_, get_objid_ext(vmg_ val2.val.obj));
6048
6049 /* discard the GC protection items */
6050 G_stk->discard(2);
6051 }
6052
6053 /*
6054 * negate the number
6055 */
neg_val(VMG_ vm_val_t * result,vm_obj_id_t self)6056 void CVmObjBigNum::neg_val(VMG_ vm_val_t *result, vm_obj_id_t self)
6057 {
6058 char *new_ext;
6059 size_t prec = get_prec(ext_);
6060
6061 /*
6062 * If I'm not an ordinary number or an infinity, return myself
6063 * unchanged. Note that we change sign for an infinity, even though
6064 * this might not make a great deal of sense mathematically.
6065 *
6066 * If I'm zero, likewise return myself unchanged. Negative zero is
6067 * still zero.
6068 */
6069 if ((get_type(ext_) != VMBN_T_NUM && get_type(ext_) != VMBN_T_INF)
6070 || is_zero(ext_))
6071 {
6072 /* return myself unchanged */
6073 result->set_obj(self);
6074 return;
6075 }
6076
6077 /* push a self-reference while we're working */
6078 G_stk->push()->set_obj(self);
6079
6080 /* if I'm an infinity, we don't need any precision in the result */
6081 if (get_type(ext_) == VMBN_T_INF)
6082 prec = 1;
6083
6084 /* create a new number with the same precision as the original */
6085 result->set_obj(create(vmg_ FALSE, prec));
6086 new_ext = get_objid_ext(vmg_ result->val.obj);
6087
6088 /* make a copy in the new object */
6089 memcpy(new_ext, ext_, calc_alloc(prec));
6090
6091 /* reverse the sign */
6092 set_neg(new_ext, !get_neg(new_ext));
6093
6094 /* remove my self-reference */
6095 G_stk->discard();
6096 }
6097
6098 /* ------------------------------------------------------------------------ */
6099 /*
6100 * check a value for equality
6101 */
equals(VMG_ vm_obj_id_t self,const vm_val_t * val,int) const6102 int CVmObjBigNum::equals(VMG_ vm_obj_id_t self, const vm_val_t *val,
6103 int /*depth*/) const
6104 {
6105 vm_val_t val2;
6106 int ret;
6107
6108 /* if the other value is a reference to self, we certainly match */
6109 if (val->typ == VM_OBJ && val->val.obj == self)
6110 return TRUE;
6111
6112 /* convert it */
6113 val2 = *val;
6114 if (!cvt_to_bignum(vmg_ self, &val2))
6115 {
6116 /* this type is not convertible to BigNumber - it's not equal */
6117 return FALSE;
6118 }
6119
6120 /* push our values for safekeeping from the garbage collector */
6121 G_stk->push()->set_obj(self);
6122 G_stk->push(&val2);
6123
6124 /* check for equality and return the result */
6125 ret = compute_eq_exact(ext_, get_objid_ext(vmg_ val2.val.obj));
6126
6127 /* discard the stacked values */
6128 G_stk->discard(2);
6129
6130 /* return the result */
6131 return ret;
6132 }
6133
6134 /*
6135 * Hash value calculation
6136 */
calc_hash(VMG_ vm_obj_id_t self,int) const6137 uint CVmObjBigNum::calc_hash(VMG_ vm_obj_id_t self, int /*depth*/) const
6138 {
6139 uint i;
6140 uint hash;
6141
6142 /* add up the digits in the number */
6143 for (hash = 0, i = 0 ; i < get_prec(ext_) ; ++i)
6144 {
6145 /* add this digit into the hash so far */
6146 hash += get_dig(ext_, i);
6147 }
6148
6149 /* add in the exponent as well */
6150 hash += (uint)get_exp(ext_);
6151
6152 /* return the combined hash */
6153 return hash;
6154 }
6155
6156 /*
6157 * compare to another value
6158 */
compare_to(VMG_ vm_obj_id_t self,const vm_val_t * val) const6159 int CVmObjBigNum::compare_to(VMG_ vm_obj_id_t self,
6160 const vm_val_t *val) const
6161 {
6162 vm_val_t val2;
6163 const char *ext2;
6164 int ret;
6165
6166 /* convert it */
6167 val2 = *val;
6168 if (!cvt_to_bignum(vmg_ self, &val2))
6169 {
6170 /* this type is not convertible to BigNumber - it's not comparable */
6171 err_throw(VMERR_INVALID_COMPARISON);
6172 }
6173
6174 /* get the other object's extension */
6175 ext2 = get_objid_ext(vmg_ val2.val.obj);
6176
6177 /* if either one is not a number, they can't be compared */
6178 if (is_nan(ext_) || is_nan(ext2))
6179 err_throw(VMERR_INVALID_COMPARISON);
6180
6181 /* if the signs differ, the positive one is greater */
6182 if (get_neg(ext_) != get_neg(ext2))
6183 {
6184 /*
6185 * if I'm negative, I'm the lesser number; otherwise I'm the
6186 * greater number
6187 */
6188 return (get_neg(ext_) ? -1 : 1);
6189 }
6190
6191 /* the signs are the same - compare the absolute values */
6192 ret = compare_abs(ext_, ext2);
6193
6194 /*
6195 * If the numbers are negative, and my absolute value is greater,
6196 * I'm actually the lesser number; if they're negative and my
6197 * absolute value is lesser, I'm the greater number. So, if I'm
6198 * negative, invert the sense of the result. Otherwise, both
6199 * numbers are positive, so the sense of the absolute value
6200 * comparison is the same as the sense of the actual comparison, so
6201 * just return the result directly.
6202 */
6203 return (get_neg(ext_) ? -ret : ret);
6204 }
6205
6206
6207 /* ------------------------------------------------------------------------ */
6208 /*
6209 * Initialize for a computation involving two operands. Checks the
6210 * operands for non-number values; if either is NAN or INF, allocates a
6211 * result value that is the same non-number type and returns null. If
6212 * both are valid numbers, we'll allocate a result value with precision
6213 * equal to the greater of the precisions of the operands, and we'll
6214 * return a pointer to the new object's extension buffer.
6215 */
compute_init_2op(VMG_ vm_val_t * result,const char * ext1,const char * ext2)6216 char *CVmObjBigNum::compute_init_2op(VMG_ vm_val_t *result,
6217 const char *ext1, const char *ext2)
6218 {
6219 size_t new_prec;
6220 char *new_ext;
6221
6222 /* get the greater precision - this is the precision of the result */
6223 new_prec = get_prec(ext1);
6224 if (get_prec(ext2) > new_prec)
6225 new_prec = get_prec(ext2);
6226
6227 /*
6228 * if either operand is not an ordinary number, we need minimal
6229 * precision to represent the result, since we don't actually need
6230 * to store a number
6231 */
6232 if (is_nan(ext1) || is_nan(ext2))
6233 new_prec = 1;
6234
6235 /* allocate a new object with the required precision */
6236 result->set_obj(create(vmg_ FALSE, new_prec));
6237
6238 /* get the extension buffer */
6239 new_ext = get_objid_ext(vmg_ result->val.obj);
6240
6241 /* check the first value for NAN/INF conditions */
6242 if (get_type(ext1) != VMBN_T_NUM)
6243 {
6244 /* set the result to the same non-number type and sign */
6245 set_type(new_ext, get_type(ext1));
6246 set_neg(new_ext, get_neg(ext1));
6247
6248 /* indicate that no further calculation is necessary */
6249 return 0;
6250 }
6251
6252 /* check the second number for NAN/INF conditions */
6253 if (get_type(ext2) != VMBN_T_NUM)
6254 {
6255 /* set the result to the same non-number type and sign */
6256 set_type(new_ext, get_type(ext2));
6257 set_neg(new_ext, get_neg(ext2));
6258
6259 /* indicate that no further calculation is necessary */
6260 return 0;
6261 }
6262
6263 /* the operands are valid - return the new extension buffer */
6264 return new_ext;
6265 }
6266
6267 /*
6268 * Compute a sum
6269 */
compute_sum(VMG_ vm_val_t * result,const char * ext1,const char * ext2)6270 void CVmObjBigNum::compute_sum(VMG_ vm_val_t *result,
6271 const char *ext1, const char *ext2)
6272 {
6273 char *new_ext;
6274
6275 /* allocate our result value */
6276 new_ext = compute_init_2op(vmg_ result, ext1, ext2);
6277
6278 /* we're done if we had a non-number operand */
6279 if (new_ext == 0)
6280 return;
6281
6282 /* compute the sum into the result */
6283 compute_sum_into(new_ext, ext1, ext2);
6284 }
6285
6286 /*
6287 * Compute a sum into the given buffer
6288 */
compute_sum_into(char * new_ext,const char * ext1,const char * ext2)6289 void CVmObjBigNum::compute_sum_into(char *new_ext,
6290 const char *ext1, const char *ext2)
6291 {
6292 /* check to see if the numbers have the same sign */
6293 if (get_neg(ext1) == get_neg(ext2))
6294 {
6295 /*
6296 * the two numbers have the same sign, so the sum has the same
6297 * sign as the two values, and the magnitude is the sum of the
6298 * absolute values of the operands
6299 */
6300
6301 /* compute the sum of the absolute values */
6302 compute_abs_sum_into(new_ext, ext1, ext2);
6303
6304 /* set the sign to match that of the operand */
6305 set_neg(new_ext, get_neg(ext1));
6306 }
6307 else
6308 {
6309 /*
6310 * one is positive and the other is negative - subtract the
6311 * absolute value of the negative one from the absolute value of
6312 * the positive one; the sign will be set correctly by the
6313 * differencing
6314 */
6315 if (get_neg(ext1))
6316 compute_abs_diff_into(new_ext, ext2, ext1);
6317 else
6318 compute_abs_diff_into(new_ext, ext1, ext2);
6319 }
6320 }
6321
6322 /*
6323 * Compute a difference
6324 */
compute_diff(VMG_ vm_val_t * result,const char * ext1,const char * ext2)6325 void CVmObjBigNum::compute_diff(VMG_ vm_val_t *result,
6326 const char *ext1, const char *ext2)
6327 {
6328 char *new_ext;
6329
6330 /* allocate our result value */
6331 new_ext = compute_init_2op(vmg_ result, ext1, ext2);
6332
6333 /* we're done if we had a non-number operand */
6334 if (new_ext == 0)
6335 return;
6336
6337 /* check to see if the numbers have the same sign */
6338 if (get_neg(ext1) == get_neg(ext2))
6339 {
6340 /*
6341 * The two numbers have the same sign, so the difference is the
6342 * difference in the magnitudes, and has a sign depending on the
6343 * order of the difference and the signs of the original values.
6344 * If both values are positive, the difference is negative if
6345 * the second value is larger than the first. If both values
6346 * are negative, the difference is negative if the second value
6347 * has larger absolute value than the first.
6348 */
6349
6350 /* compute the difference in magnitudes */
6351 compute_abs_diff_into(new_ext, ext1, ext2);
6352
6353 /*
6354 * if the original values were negative, then the sign of the
6355 * result is the opposite of the sign of the difference of the
6356 * absolute values; otherwise, it's the same as the sign of the
6357 * difference of the absolute values
6358 */
6359 if (get_neg(ext1))
6360 negate(new_ext);
6361 }
6362 else
6363 {
6364 /*
6365 * one is positive and the other is negative - the result has
6366 * the sign of the first operand, and has magnitude equal to the
6367 * sum of the absolute values
6368 */
6369
6370 /* compute the sum of the absolute values */
6371 compute_abs_sum_into(new_ext, ext1, ext2);
6372
6373 /* set the sign of the result to that of the first operand */
6374 if (!is_zero(new_ext))
6375 set_neg(new_ext, get_neg(ext1));
6376 }
6377 }
6378
6379 /*
6380 * Compute a product
6381 */
compute_prod(VMG_ vm_val_t * result,const char * ext1,const char * ext2)6382 void CVmObjBigNum::compute_prod(VMG_ vm_val_t *result,
6383 const char *ext1, const char *ext2)
6384 {
6385 char *new_ext;
6386
6387 /* allocate our result value */
6388 new_ext = compute_init_2op(vmg_ result, ext1, ext2);
6389
6390 /* we're done if we had a non-number operand */
6391 if (new_ext == 0)
6392 return;
6393
6394 /* compute the product */
6395 compute_prod_into(new_ext, ext1, ext2);
6396 }
6397
6398 /*
6399 * Compute a quotient
6400 */
compute_quotient(VMG_ vm_val_t * result,const char * ext1,const char * ext2)6401 void CVmObjBigNum::compute_quotient(VMG_ vm_val_t *result,
6402 const char *ext1, const char *ext2)
6403 {
6404 char *new_ext;
6405
6406 /* allocate our result value */
6407 new_ext = compute_init_2op(vmg_ result, ext1, ext2);
6408
6409 /* we're done if we had a non-number operand */
6410 if (new_ext == 0)
6411 return;
6412
6413 /* compute the quotient */
6414 compute_quotient_into(vmg_ new_ext, 0, ext1, ext2);
6415 }
6416
6417 /*
6418 * Determine if two values are equal, rounding the value with greater
6419 * precision to the shorter precision if the two are not of equal
6420 * precision
6421 */
compute_eq_round(VMG_ const char * ext1,const char * ext2)6422 int CVmObjBigNum::compute_eq_round(VMG_ const char *ext1, const char *ext2)
6423 {
6424 size_t prec1 = get_prec(ext1);
6425 size_t prec2 = get_prec(ext2);
6426 const char *shorter;
6427 const char *longer;
6428 char *tmp_ext;
6429 uint tmp_hdl;
6430 int ret;
6431
6432 /*
6433 * allocate a temporary register with a rounded copy of the more
6434 * precise of the values
6435 */
6436 if (prec1 > prec2)
6437 {
6438 /* the first one is longer */
6439 longer = ext1;
6440 shorter = ext2;
6441 }
6442 else if (prec1 < prec2)
6443 {
6444 /* the second one is longer */
6445 longer = ext2;
6446 shorter = ext1;
6447 }
6448 else
6449 {
6450 /* they're the same - do an exact comparison */
6451 return compute_eq_exact(ext1, ext2);
6452 }
6453
6454 /* get a temp register for rounding the longer value */
6455 alloc_temp_regs(vmg_ get_prec(shorter), 1, &tmp_ext, &tmp_hdl);
6456
6457 /* make a rounded copy */
6458 copy_val(tmp_ext, longer, TRUE);
6459
6460 /* compare the rounded copy of the longer value to the shorter value */
6461 ret = compute_eq_exact(shorter, tmp_ext);
6462
6463 /* release the temporary register */
6464 release_temp_regs(vmg_ 1, tmp_hdl);
6465
6466 /* return the result */
6467 return ret;
6468 }
6469
6470 /*
6471 * Make an exact comparison for equality. If one value is more precise
6472 * than the other, we'll implicitly extend the shorter value to the
6473 * right with trailing zeroes.
6474 */
compute_eq_exact(const char * ext1,const char * ext2)6475 int CVmObjBigNum::compute_eq_exact(const char *ext1, const char *ext2)
6476 {
6477 const char *longer;
6478 size_t min_prec;
6479 size_t max_prec;
6480 size_t prec1;
6481 size_t prec2;
6482 size_t idx;
6483
6484 /*
6485 * if either is not an ordinary number, they are never equal to any
6486 * other value (note that this means INF != INF and NAN != NAN,
6487 * which is reasonable because these values cannot be meaningfully
6488 * compared; one NAN might mean something totally different from
6489 * another, and likewise various infinities are not comparable)
6490 */
6491 if (is_nan(ext1) || is_nan(ext2))
6492 return FALSE;
6493
6494 /* figure out if one is more precise than the other */
6495 prec1 = get_prec(ext1);
6496 prec2 = get_prec(ext2);
6497 if (prec1 > prec2)
6498 {
6499 /* ext1 is longer */
6500 longer = ext1;
6501 max_prec = prec1;
6502 min_prec = prec2;
6503 }
6504 else if (prec2 > prec1)
6505 {
6506 /* ext2 is longer */
6507 longer = ext2;
6508 max_prec = prec2;
6509 min_prec = prec1;
6510 }
6511 else
6512 {
6513 /* they're the same */
6514 longer = 0;
6515 min_prec = max_prec = prec1;
6516 }
6517
6518 /* if the signs aren't the same, the numbers are not equal */
6519 if (get_neg(ext1) != get_neg(ext2))
6520 return FALSE;
6521
6522 /* if the exponents aren't equal, the numbers are not equal */
6523 if (get_exp(ext1) != get_exp(ext2))
6524 return FALSE;
6525
6526 /*
6527 * compare digits up to the smaller precision, then make sure that
6528 * the larger-precision value's digits are all zeroes from there out
6529 */
6530 for (idx = 0 ; idx < min_prec ; ++idx)
6531 {
6532 /* if they don't match, return not-equal */
6533 if (get_dig(ext1, idx) != get_dig(ext2, idx))
6534 return FALSE;
6535 }
6536
6537 /* check the longer one to make sure it's all zeroes */
6538 if (longer != 0)
6539 {
6540 /* scan the remainder of the longer one */
6541 for ( ; idx < max_prec ; ++idx)
6542 {
6543 /* if this digit is non-zero, it's not a match */
6544 if (get_dig(longer, idx) != 0)
6545 return FALSE;
6546 }
6547 }
6548
6549 /* it's a match */
6550 return TRUE;
6551 }
6552
6553 /* ------------------------------------------------------------------------ */
6554 /*
6555 * Compute the sum of two absolute values into the given buffer
6556 */
compute_abs_sum_into(char * new_ext,const char * ext1,const char * ext2)6557 void CVmObjBigNum::compute_abs_sum_into(char *new_ext,
6558 const char *ext1, const char *ext2)
6559 {
6560 int max_exp;
6561 int lo1, hi1;
6562 int lo2, hi2;
6563 int lo3, hi3;
6564 int pos;
6565 int carry;
6566 int trail_dig;
6567
6568 /* if one or the other is identically zero, return the other value */
6569 if (is_zero(ext1))
6570 {
6571 /* ext1 is zero - return ext2 */
6572 copy_val(new_ext, ext2, TRUE);
6573 return;
6574 }
6575 else if (is_zero(ext2))
6576 {
6577 /* ext2 is zero - return ext1 */
6578 copy_val(new_ext, ext1, TRUE);
6579 return;
6580 }
6581
6582 /*
6583 * start the new value with the larger of the two exponents - this
6584 * will have the desired effect of dropping the least significant
6585 * digits if any digits must be dropped
6586 */
6587 max_exp = get_exp(ext1);
6588 if (get_exp(ext2) > max_exp)
6589 max_exp = get_exp(ext2);
6590 set_exp(new_ext, max_exp);
6591
6592 /* compute the digit positions at the bounds of each of our values */
6593 hi1 = get_exp(ext1) - 1;
6594 lo1 = get_exp(ext1) - get_prec(ext1);
6595
6596 hi2 = get_exp(ext2) - 1;
6597 lo2 = get_exp(ext2) - get_prec(ext2);
6598
6599 hi3 = get_exp(new_ext) - 1;
6600 lo3 = get_exp(new_ext) - get_prec(new_ext);
6601
6602 /*
6603 * If one of the values provides a digit one past the end of our
6604 * result, remember that as the trailing digit that we're going to
6605 * drop. We'll check this when we're done to see if we need to
6606 * round the number. Since the result has precision at least as
6607 * large as the larger of the two inputs, we can only be dropping
6608 * significant digits from one of the two inputs - we can't be
6609 * cutting off both inputs.
6610 */
6611 trail_dig = 0;
6612 if (lo3-1 >= lo1 && lo3-1 <= hi1)
6613 {
6614 /* remember the digit */
6615 trail_dig = get_dig(ext1, get_exp(ext1) - (lo3-1) - 1);
6616 }
6617 else if (lo3-1 >= lo2 && lo3-1 <= hi2)
6618 {
6619 /* remember the digit */
6620 trail_dig = get_dig(ext2, get_exp(ext2) - (lo3-1) - 1);
6621 }
6622
6623 /* no carry yet */
6624 carry = 0;
6625
6626 /* add the digits */
6627 for (pos = lo3 ; pos <= hi3 ; ++pos)
6628 {
6629 int acc;
6630
6631 /* start with the carry */
6632 acc = carry;
6633
6634 /* add the first value digit if it's in range */
6635 if (pos >= lo1 && pos <= hi1)
6636 acc += get_dig(ext1, get_exp(ext1) - pos - 1);
6637
6638 /* add the second value digit if it's in range */
6639 if (pos >= lo2 && pos <= hi2)
6640 acc += get_dig(ext2, get_exp(ext2) - pos - 1);
6641
6642 /* check for carry */
6643 if (acc > 9)
6644 {
6645 /* reduce the accumulator and set the carry */
6646 acc -= 10;
6647 carry = 1;
6648 }
6649 else
6650 {
6651 /* no carry */
6652 carry = 0;
6653 }
6654
6655 /* set the digit in the result */
6656 set_dig(new_ext, get_exp(new_ext) - pos - 1, acc);
6657 }
6658
6659 /*
6660 * If we have a carry at the end, we must carry it out to a new
6661 * digit. In order to do this, we must shift the whole number right
6662 * one place, then insert the one.
6663 */
6664 if (carry)
6665 {
6666 /*
6667 * remember the last digit of the result, which we won't have
6668 * space to store after the shift
6669 */
6670 trail_dig = get_dig(new_ext, get_prec(new_ext) - 1);
6671
6672 /* shift the result right */
6673 shift_right(new_ext, 1);
6674
6675 /* increase the exponent to compensate for the shift */
6676 set_exp(new_ext, get_exp(new_ext) + 1);
6677
6678 /* set the leading 1 */
6679 set_dig(new_ext, 0, 1);
6680 }
6681
6682 /* the sum of two absolute values is always positive */
6683 set_neg(new_ext, FALSE);
6684
6685 /* round up the value if the trailing digit is 5 or higher */
6686 if (trail_dig >= 5)
6687 round_up_abs(new_ext);
6688
6689 /* normalize the number */
6690 normalize(new_ext);
6691 }
6692
6693 /*
6694 * Compute the difference of two absolute values into the given buffer
6695 */
compute_abs_diff_into(char * new_ext,const char * ext1,const char * ext2)6696 void CVmObjBigNum::compute_abs_diff_into(char *new_ext,
6697 const char *ext1, const char *ext2)
6698 {
6699 int max_exp;
6700 int lo1, hi1;
6701 int lo2, hi2;
6702 int lo3, hi3;
6703 int pos;
6704 int result_neg = FALSE;
6705 int borrow;
6706
6707 /* if we're subtracting zero or from zero, it's easy */
6708 if (is_zero(ext2))
6709 {
6710 /*
6711 * we're subtracting zero from another value - the result is
6712 * simply the first value
6713 */
6714 copy_val(new_ext, ext1, TRUE);
6715 return;
6716 }
6717 else if (is_zero(ext1))
6718 {
6719 /*
6720 * We're subtracting a value from zero - we know the value we're
6721 * subtracting is non-zero (we already checked for that case
6722 * above), and we're only considering the absolute values, so
6723 * simply return the negative of the absolute value of the
6724 * second operand.
6725 */
6726 copy_val(new_ext, ext2, TRUE);
6727 set_neg(new_ext, TRUE);
6728 return;
6729 }
6730
6731 /*
6732 * Compare the absolute values of the two operands. If the first
6733 * value is larger than the second, subtract them in the given
6734 * order. Otherwise, reverse the order and note that the result is
6735 * negative.
6736 */
6737 if (compare_abs(ext1, ext2) < 0)
6738 {
6739 const char *tmp;
6740
6741 /* the result will be negative */
6742 result_neg = TRUE;
6743
6744 /* swap the order of the subtraction */
6745 tmp = ext1;
6746 ext1 = ext2;
6747 ext2 = tmp;
6748 }
6749
6750 /*
6751 * start the new value with the larger of the two exponents - this
6752 * will have the desired effect of dropping the least significant
6753 * digits if any digits must be dropped
6754 */
6755 max_exp = get_exp(ext1);
6756 if (get_exp(ext2) > max_exp)
6757 max_exp = get_exp(ext2);
6758 set_exp(new_ext, max_exp);
6759
6760 /* compute the digit positions at the bounds of each of our values */
6761 hi1 = get_exp(ext1) - 1;
6762 lo1 = get_exp(ext1) - get_prec(ext1);
6763
6764 hi2 = get_exp(ext2) - 1;
6765 lo2 = get_exp(ext2) - get_prec(ext2);
6766
6767 hi3 = get_exp(new_ext) - 1;
6768 lo3 = get_exp(new_ext) - get_prec(new_ext);
6769
6770 /* start off with no borrow */
6771 borrow = 0;
6772
6773 /*
6774 * Check for borrowing from before the least significant digit
6775 * position in common to both numbers
6776 */
6777 if (lo3-1 >= lo1 && lo3-1 <= hi1)
6778 {
6779 /*
6780 * In this case, we would be dropping precision from the end of
6781 * the top number. This case should never happen - the only way
6782 * it could happen is for the bottom number to extend to the
6783 * left of the top number at the most significant end. This
6784 * cannot happen because the bottom number always has small
6785 * magnitude by this point (we guarantee this above). So, we
6786 * should never get here.
6787 */
6788 assert(FALSE);
6789 }
6790 else if (lo3-1 >= lo2 && lo3-1 <= hi2)
6791 {
6792 size_t idx;
6793
6794 /*
6795 * We're dropping precision from the bottom number, so we want
6796 * to borrow into the subtraction if the rest of the number is
6797 * greater than 5xxxx. If it's exactly 5000..., do not borrow,
6798 * since the result would end in 5 and we'd round up.
6799 *
6800 * If the next digit is 6 or greater, we know for a fact we'll
6801 * have to borrow. If the next digit is 4 or less, we know for
6802 * a fact we won't have to borrow. If the next digit is 5,
6803 * though, we must look at the rest of the number to see if
6804 * there's anything but trailing zeroes.
6805 */
6806 idx = (size_t)(get_exp(ext2) - (lo3-1) - 1);
6807 if (get_dig(ext2, idx) >= 6)
6808 {
6809 /* definitely borrow */
6810 borrow = 1;
6811 }
6812 else if (get_dig(ext2, idx) == 5)
6813 {
6814 /* borrow only if we have something non-zero following */
6815 for (++idx ; idx < get_prec(ext2) ; ++idx)
6816 {
6817 /* if it's non-zero, we must borrow */
6818 if (get_dig(ext2, idx) != 0)
6819 {
6820 /* note the borrow */
6821 borrow = 1;
6822
6823 /* no need to keep scanning */
6824 break;
6825 }
6826 }
6827 }
6828 }
6829
6830 /* subtract the digits from least to most significant */
6831 for (pos = lo3 ; pos <= hi3 ; ++pos)
6832 {
6833 int acc;
6834
6835 /* start with zero in the accumulator */
6836 acc = 0;
6837
6838 /* start with the top-line value if it's represented here */
6839 if (pos >= lo1 && pos <= hi1)
6840 acc = get_dig(ext1, get_exp(ext1) - pos - 1);
6841
6842 /* subtract the second value digit if it's represented here */
6843 if (pos >= lo2 && pos <= hi2)
6844 acc -= get_dig(ext2, get_exp(ext2) - pos - 1);
6845
6846 /* subtract the borrow */
6847 acc -= borrow;
6848
6849 /* check for borrow */
6850 if (acc < 0)
6851 {
6852 /* increase the accumulator */
6853 acc += 10;
6854
6855 /* we must borrow from the next digit up */
6856 borrow = 1;
6857 }
6858 else
6859 {
6860 /* we're in range - no need to borrow */
6861 borrow = 0;
6862 }
6863
6864 /* set the digit in the result */
6865 set_dig(new_ext, get_exp(new_ext) - pos - 1, acc);
6866 }
6867
6868 /* set the sign of the result as calculated earlier */
6869 set_neg(new_ext, result_neg);
6870
6871 /* normalize the number */
6872 normalize(new_ext);
6873 }
6874
6875 /*
6876 * Compute the product of the values into the given buffer
6877 */
compute_prod_into(char * new_ext,const char * ext1,const char * ext2)6878 void CVmObjBigNum::compute_prod_into(char *new_ext,
6879 const char *ext1, const char *ext2)
6880 {
6881 size_t prec1 = get_prec(ext1);
6882 size_t prec2 = get_prec(ext2);
6883 size_t new_prec = get_prec(new_ext);
6884 size_t idx1;
6885 size_t idx2;
6886 size_t out_idx;
6887 size_t start_idx;
6888 int out_exp;
6889 int trail_dig;
6890
6891 /* start out with zero in the accumulator */
6892 memset(new_ext + VMBN_MANT, 0, (new_prec + 1)/2);
6893
6894 /*
6895 * Initially write output in the same 'column' where the top number
6896 * ends, so we start out with the same scale as the top number.
6897 */
6898 start_idx = get_prec(ext1);
6899
6900 /*
6901 * initial result exponent is the sum of the exponents, minus the
6902 * number of digits in the bottom number (effectively, this lets us
6903 * treat the bottom number as a whole number by scaling it enough to
6904 * make it whole, soaking up the factors of ten into decimal point
6905 * relocation)
6906 */
6907 out_exp = get_exp(ext1) + get_exp(ext2) - prec2;
6908
6909 /* there's no trailing accumulator digit yet */
6910 trail_dig = 0;
6911
6912 /*
6913 * Loop over digits in the bottom number, from least significant to
6914 * most significant - we'll multiply each digit of the bottom number
6915 * by the top number and add the result into the accumulator.
6916 */
6917 for (idx2 = prec2 ; idx2 != 0 ; )
6918 {
6919 int carry;
6920 int dig;
6921 int ext2_dig;
6922
6923 /* no carry yet on this round */
6924 carry = 0;
6925
6926 /* start writing this round at the output start index */
6927 out_idx = start_idx;
6928
6929 /* move to the next digit */
6930 --idx2;
6931
6932 /* get this digit of ext2 */
6933 ext2_dig = get_dig(ext2, idx2);
6934
6935 /*
6936 * if this digit of ext2 is non-zero, multiply the digits of
6937 * ext1 by the digit (obviously if the digit is zero, there's no
6938 * need to iterate over the digits of ext1)
6939 */
6940 if (ext2_dig != 0)
6941 {
6942 /*
6943 * loop over digits in the top number, from least to most
6944 * significant
6945 */
6946 for (idx1 = prec1 ; idx1 != 0 ; )
6947 {
6948 /* move to the next digit of the top number */
6949 --idx1;
6950
6951 /* move to the next digit of the accumulator */
6952 --out_idx;
6953
6954 /*
6955 * compute the product of the current digits, and add in
6956 * the carry from the last digit, then add in the
6957 * current accumulator digit in this position
6958 */
6959 dig = (get_dig(ext1, idx1) * ext2_dig)
6960 + carry
6961 + get_dig(new_ext, out_idx);
6962
6963 /* figure the carry to the next digit */
6964 carry = (dig / 10);
6965 dig = dig % 10;
6966
6967 /* store the new digit */
6968 set_dig(new_ext, out_idx, dig);
6969 }
6970 }
6971
6972 /*
6973 * Shift the result accumulator right in preparation for the
6974 * next round. One exception: if this is the last (most
6975 * significant) digit of the bottom number, and there's no
6976 * carry, there's no need to shift the number, since we'd just
6977 * normalize away the leading zero anyway
6978 */
6979 if (idx2 != 0 || carry != 0)
6980 {
6981 /* remember the trailing digit that we're going to drop */
6982 trail_dig = get_dig(new_ext, new_prec - 1);
6983
6984 /* shift the accumulator */
6985 shift_right(new_ext, 1);
6986
6987 /* increase the output exponent */
6988 ++out_exp;
6989
6990 /* insert the carry as the lead digit */
6991 set_dig(new_ext, 0, carry);
6992 }
6993 }
6994
6995 /* set the result exponent */
6996 set_exp(new_ext, out_exp);
6997
6998 /*
6999 * set the sign - if both inputs have the same sign, the output is
7000 * positive, otherwise it's negative
7001 */
7002 set_neg(new_ext, get_neg(ext1) != get_neg(ext2));
7003
7004 /* if the trailing digit is 5 or greater, round up */
7005 if (trail_dig >= 5)
7006 round_up_abs(new_ext);
7007
7008 /* normalize the number */
7009 normalize(new_ext);
7010 }
7011
7012 /*
7013 * Compute a quotient into the given buffer. If new_rem_ext is
7014 * non-null, we'll save the remainder into this buffer. We calculate
7015 * the remainder only to the precision of the quotient.
7016 */
compute_quotient_into(VMG_ char * new_ext,char * new_rem_ext,const char * ext1,const char * ext2)7017 void CVmObjBigNum::compute_quotient_into(VMG_ char *new_ext,
7018 char *new_rem_ext,
7019 const char *ext1, const char *ext2)
7020 {
7021 char *rem_ext;
7022 uint rem_hdl;
7023 char *rem_ext2;
7024 uint rem_hdl2;
7025 int quo_exp;
7026 size_t quo_idx;
7027 size_t quo_prec = get_prec(new_ext);
7028 size_t dvd_prec = get_prec(ext1);
7029 size_t dvs_prec = get_prec(ext2);
7030 char *dvs_ext;
7031 uint dvs_hdl;
7032 char *dvs_ext2;
7033 uint dvs_hdl2;
7034 int lead_dig_set;
7035 int zero_remainder;
7036 int acc;
7037 size_t rem_prec;
7038
7039 /* if the divisor is zero, throw an error */
7040 if (is_zero(ext2))
7041 err_throw(VMERR_DIVIDE_BY_ZERO);
7042
7043 /* if the dividend is zero, the result is zero */
7044 if (is_zero(ext1))
7045 {
7046 /* set the result to zero */
7047 set_zero(new_ext);
7048
7049 /* if they want the remainder, it's zero also */
7050 if (new_rem_ext != 0)
7051 set_zero(new_rem_ext);
7052
7053 /* we're done */
7054 return;
7055 }
7056
7057 /*
7058 * Calculate the precision we need for the running remainder. We
7059 * must retain in the remainder enough precision to calculate exact
7060 * differences, so we need the greater of the precisions of the
7061 * dividend and the divisor, plus enough extra digits for the
7062 * maximum relative shifting. We will have to shift at most one
7063 * extra digit, but use two to be extra safe.
7064 */
7065 rem_prec = dvd_prec;
7066 if (rem_prec < dvs_prec)
7067 rem_prec = dvs_prec;
7068 rem_prec += 2;
7069
7070 /*
7071 * Make sure the precision is at least three, since it simplifies
7072 * our digit approximation calculation.
7073 */
7074 if (rem_prec < 3)
7075 rem_prec = 3;
7076
7077 /*
7078 * Allocate two temporary registers for the running remainder, and
7079 * one more for the multiplied divisor. Note that we allocate the
7080 * multiplied divisor at our higher precision so that we don't lose
7081 * digits in our multiplier calculations.
7082 */
7083 alloc_temp_regs(vmg_ rem_prec, 3,
7084 &rem_ext, &rem_hdl, &rem_ext2, &rem_hdl2,
7085 &dvs_ext2, &dvs_hdl2);
7086
7087 /*
7088 * Allocate another temp register for the shifted divisor. Note
7089 * that we need a different precision here, which is why we must
7090 * make a separate allocation call
7091 */
7092 err_try
7093 {
7094 /* make the additional allocation */
7095 alloc_temp_regs(vmg_ dvs_prec, 1, &dvs_ext, &dvs_hdl);
7096 }
7097 err_catch(exc)
7098 {
7099 /* delete the first group of registers we allocated */
7100 release_temp_regs(vmg_ 2, rem_hdl, rem_hdl2);
7101
7102 /* re-throw the exception */
7103 err_rethrow();
7104 }
7105 err_end;
7106
7107 /* the dividend is the initial value of the running remainder */
7108 copy_val(rem_ext, ext1, TRUE);
7109
7110 /* copy the initial divisor into our temp register */
7111 copy_val(dvs_ext, ext2, TRUE);
7112
7113 /* we haven't set a non-zero leading digit yet */
7114 lead_dig_set = FALSE;
7115
7116 /*
7117 * scale the divisor so that the divisor and dividend have the same
7118 * exponent, and absorb the multiplier in the quotient
7119 */
7120 quo_exp = get_exp(ext1) - get_exp(ext2) + 1;
7121 set_exp(dvs_ext, get_exp(ext1));
7122
7123 /* we don't have a zero remainder yet */
7124 zero_remainder = FALSE;
7125
7126 /*
7127 * if the quotient is going to be entirely fractional, the dividend
7128 * is the remainder, and the quotient is zero
7129 */
7130 if (new_rem_ext != 0 && quo_exp <= 0)
7131 {
7132 /* copy the initial remainder into the output remainder */
7133 copy_val(new_rem_ext, rem_ext, TRUE);
7134
7135 /* set the quotient to zero */
7136 set_zero(new_ext);
7137
7138 /* we have the result - no need to do any more work */
7139 goto done;
7140 }
7141
7142 /*
7143 * Loop over each digit of precision of the quotient.
7144 */
7145 for (quo_idx = 0 ; ; )
7146 {
7147 int rem_approx, dvs_approx;
7148 int dig_approx;
7149 char *tmp;
7150 int exp_diff;
7151
7152 /* start out with 0 in our digit accumulator */
7153 acc = 0;
7154
7155 /*
7156 * Get the first few digits of the remainder, and the first few
7157 * digits of the divisor, rounding up the divisor and rounding
7158 * down the remainder. Compute the quotient - this will give us
7159 * a rough guess and a lower bound for the current digit's
7160 * value.
7161 */
7162 rem_approx = (get_dig(rem_ext, 0)*100
7163 + get_dig(rem_ext, 1)*10
7164 + get_dig(rem_ext, 2));
7165 dvs_approx = (get_dig(dvs_ext, 0)*100
7166 + (dvs_prec >= 2 ? get_dig(dvs_ext, 1) : 0)*10
7167 + (dvs_prec >= 3 ? get_dig(dvs_ext, 2) : 0)
7168 + 1);
7169
7170 /* adjust for differences in the scale */
7171 exp_diff = get_exp(rem_ext) - get_exp(dvs_ext);
7172 if (exp_diff > 0)
7173 {
7174 /* the remainder is larger - scale it up */
7175 for ( ; exp_diff > 0 ; --exp_diff)
7176 rem_approx *= 10;
7177 }
7178 else if (exp_diff <= -3)
7179 {
7180 /*
7181 * The divisor is larger by more than 10^3, which means that
7182 * the result of our integer division is definitely going to
7183 * be zero, so there's no point in actually doing the
7184 * calculation. This is only a special case because, for
7185 * sufficiently large negative differences, we'd have to
7186 * multiply our divisor approximation by 10 so many times
7187 * that we'd overflow a native int type, at which point we'd
7188 * get 0 in the divisor, which would result in a
7189 * divide-by-zero. To avoid this, just put 1000 in our
7190 * divisor, which is definitely larger than anything we can
7191 * have in rem_ext at this point (since it was just three
7192 * decimal digits, after all).
7193 */
7194 dvs_approx = 1000;
7195 }
7196 else if (exp_diff < 0)
7197 {
7198 /* the divisor is larger - scale it up */
7199 for ( ; exp_diff < 0 ; ++exp_diff)
7200 dvs_approx *= 10;
7201 }
7202
7203 /* calculate our initial guess for this digit */
7204 dig_approx = rem_approx / dvs_approx;
7205
7206 /*
7207 * If this digit is above 2, it'll save us a lot of work to
7208 * subtract digit*divisor once, rather than iteratively
7209 * deducting the divisor one time per iteration. (It costs us a
7210 * little to calculate the digit*divisor product, so we don't
7211 * want to do this for very small digit values.)
7212 */
7213 if (dig_approx > 2)
7214 {
7215 /* put the approximate digit in the accumulator */
7216 acc = dig_approx;
7217
7218 /* make a copy of the divisor */
7219 copy_val(dvs_ext2, dvs_ext, FALSE);
7220
7221 /*
7222 * multiply it by the guess for the digit - we know this
7223 * will always be less than or equal to the real value
7224 * because of the way we did the rounding
7225 */
7226 mul_by_long(dvs_ext2, (ulong)dig_approx);
7227
7228 /* subtract it from the running remainder */
7229 compute_abs_diff_into(rem_ext2, rem_ext, dvs_ext2);
7230
7231 /* if that leaves zero in the remainder, note it */
7232 if (is_zero(rem_ext2))
7233 {
7234 zero_remainder = TRUE;
7235 break;
7236 }
7237
7238 /*
7239 * swap the remainder registers, since rem_ext2 is now the
7240 * new running remainder value
7241 */
7242 tmp = rem_ext;
7243 rem_ext = rem_ext2;
7244 rem_ext2 = tmp;
7245
7246 /*
7247 * Now we'll finish up the job by subtracting the divisor
7248 * from the remainder as many more times as necessary for
7249 * the remainder to fall below the divisor. We can't be
7250 * exact at this step because we're not considering all of
7251 * the digits, but we should only have one more subtraction
7252 * remaining at this point.
7253 */
7254 }
7255
7256 /*
7257 * subtract the divisor from the running remainder as many times
7258 * as we can
7259 */
7260 for ( ; ; ++acc)
7261 {
7262 int comp_res;
7263
7264 /* compare the running remainder to the divisor */
7265 comp_res = compare_abs(rem_ext, dvs_ext);
7266 if (comp_res < 0)
7267 {
7268 /*
7269 * the remainder is smaller than the divisor - we have
7270 * our result for this digit
7271 */
7272 break;
7273 }
7274 else if (comp_res == 0)
7275 {
7276 /* note that we have a zero remainder */
7277 zero_remainder = TRUE;
7278
7279 /* count one more subtraction */
7280 ++acc;
7281
7282 /* we have our result for this digit */
7283 break;
7284 }
7285
7286 /* subtract it */
7287 compute_abs_diff_into(rem_ext2, rem_ext, dvs_ext);
7288
7289 /*
7290 * swap the remainder registers, since rem_ext2 is now the
7291 * new running remainder value
7292 */
7293 tmp = rem_ext;
7294 rem_ext = rem_ext2;
7295 rem_ext2 = tmp;
7296 }
7297
7298 /* store this digit of the quotient */
7299 if (quo_idx < quo_prec)
7300 {
7301 /* store the digit */
7302 set_dig(new_ext, quo_idx, acc);
7303 }
7304 else if (quo_idx == quo_prec)
7305 {
7306 /* set the quotient's exponent */
7307 set_exp(new_ext, quo_exp);
7308
7309 /*
7310 * this is the last digit, which we calculated for rounding
7311 * purposes only - if it's 5 or greater, round up the value,
7312 * otherwise leave it as it is
7313 */
7314 if (acc >= 5)
7315 round_up_abs(new_ext);
7316
7317 /* we've reached the rounding digit - we can stop now */
7318 break;
7319 }
7320
7321 /*
7322 * if this is a non-zero digit, we've found a significant
7323 * leading digit
7324 */
7325 if (acc != 0)
7326 lead_dig_set = TRUE;
7327
7328 /*
7329 * if we've found a significant leading digit, move to the next
7330 * digit of the quotient; if not, adjust the quotient exponent
7331 * down one, and keep preparing to write the first digit at the
7332 * first "column" of the quotient
7333 */
7334 if (lead_dig_set)
7335 ++quo_idx;
7336 else
7337 --quo_exp;
7338
7339 /*
7340 * If we have an exact result (a zero remainder), we're done.
7341 *
7342 * Similarly, if we've reached the units digit, and the caller
7343 * wants the remainder, stop now - the caller wants an integer
7344 * result for the quotient, which we now have.
7345 *
7346 * Similarly, if we've reached the rounding digit, and the
7347 * caller wants the remainder, skip the rounding step - the
7348 * caller wants an unrounded result for the quotient so that the
7349 * quotient times the divisor plus the remainder equals the
7350 * dividend.
7351 */
7352 if (zero_remainder
7353 || (new_rem_ext != 0
7354 && ((int)quo_idx == quo_exp || quo_idx == quo_prec)))
7355 {
7356 /* zero any remaining digits of the quotient */
7357 for ( ; quo_idx < quo_prec ; ++quo_idx)
7358 set_dig(new_ext, quo_idx, 0);
7359
7360 /* set the quotient's exponent */
7361 set_exp(new_ext, quo_exp);
7362
7363 /* that's it */
7364 break;
7365 }
7366
7367 /* divide the divisor by ten */
7368 set_exp(dvs_ext, get_exp(dvs_ext) - 1);
7369 }
7370
7371 /* store the remainder if the caller wants the value */
7372 if (new_rem_ext != 0)
7373 {
7374 /* save the remainder into the caller's buffer */
7375 if (zero_remainder)
7376 {
7377 /* the remainder is exactly zero */
7378 set_zero(new_rem_ext);
7379 }
7380 else
7381 {
7382 /* copy the running remainder */
7383 copy_val(new_rem_ext, rem_ext, TRUE);
7384
7385 /* the remainder has the same sign as the dividend */
7386 set_neg(new_rem_ext, get_neg(ext1));
7387
7388 /* normalize the remainder */
7389 normalize(new_rem_ext);
7390 }
7391 }
7392
7393 /*
7394 * the quotient is positive if both the divisor and dividend have
7395 * the same sign, negative if they're different
7396 */
7397 set_neg(new_ext, get_neg(ext1) != get_neg(ext2));
7398
7399 /* normalize the quotient */
7400 normalize(new_ext);
7401
7402 done:
7403 /* release the temporary registers */
7404 release_temp_regs(vmg_ 4, rem_hdl, rem_hdl2, dvs_hdl, dvs_hdl2);
7405 }
7406
7407 /*
7408 * Compare the absolute values of two numbers
7409 */
compare_abs(const char * ext1,const char * ext2)7410 int CVmObjBigNum::compare_abs(const char *ext1, const char *ext2)
7411 {
7412 size_t idx;
7413 int zero1 = is_zero(ext1);
7414 int zero2 = is_zero(ext2);
7415 size_t prec1 = get_prec(ext1);
7416 size_t prec2 = get_prec(ext2);
7417 size_t max_prec;
7418 size_t min_prec;
7419 const char *max_ext;
7420 int result;
7421
7422 /*
7423 * if one is zero and the other is not, the one that's not zero has
7424 * a larger absolute value
7425 */
7426 if (zero1)
7427 return (zero2 ? 0 : -1);
7428 if (zero2)
7429 return (zero1 ? 0 : 1);
7430
7431 /*
7432 * if the exponents differ, the one with the larger exponent is the
7433 * larger number (this is necessarily true because we keep all
7434 * numbers normalized)
7435 */
7436 if (get_exp(ext1) > get_exp(ext2))
7437 return 1;
7438 else if (get_exp(ext1) < get_exp(ext2))
7439 return -1;
7440
7441 /*
7442 * The exponents are equal, so we must compare the mantissas digit
7443 * by digit
7444 */
7445
7446 /* get the larger of the two precisions */
7447 min_prec = prec2;
7448 max_prec = prec1;
7449 max_ext = ext1;
7450 if (prec2 > max_prec)
7451 {
7452 min_prec = prec1;
7453 max_prec = prec2;
7454 max_ext = ext2;
7455 }
7456
7457 /*
7458 * The digits are in order from most significant to least
7459 * significant, which means we can use memcmp to compare the common
7460 * parts. However, we can only compare an even number of digits
7461 * this way, so round down the common precision if it's odd.
7462 */
7463 if (min_prec > 1
7464 && (result = memcmp(ext1 + VMBN_MANT, ext2 + VMBN_MANT,
7465 min_prec/2)) != 0)
7466 {
7467 /*
7468 * they're different in the common memcmp-able parts, so return
7469 * the memcmp result
7470 */
7471 return result;
7472 }
7473
7474 /* if the common precision is odd, compare the final common digit */
7475 if ((min_prec & 1) != 0
7476 && (result = ((int)get_dig(ext1, min_prec-1)
7477 - (int)get_dig(ext2, min_prec-1))) != 0)
7478 return result;
7479
7480 /*
7481 * the common parts of the mantissas are identical; check each
7482 * remaining digit of the longer one to see if any are non-zero
7483 */
7484 for (idx = min_prec ; idx < max_prec ; ++idx)
7485 {
7486 /*
7487 * if this digit is non-zero, the longer one is greater, because
7488 * the shorter one has an implied zero in this position
7489 */
7490 if (get_dig(max_ext, idx) != 0)
7491 return (int)prec1 - (int)prec2;
7492 }
7493
7494 /* they're identical */
7495 return 0;
7496 }
7497
7498 /* ------------------------------------------------------------------------ */
7499 /*
7500 * Round a value
7501 */
round_val(VMG_ vm_val_t * new_val,const char * ext,size_t digits,int always_create)7502 const char *CVmObjBigNum::round_val(VMG_ vm_val_t *new_val, const char *ext,
7503 size_t digits, int always_create)
7504 {
7505 char *new_ext;
7506 int idx;
7507 int need_carry;
7508 int need_round;
7509
7510 /* presume we need rounding */
7511 need_round = TRUE;
7512
7513 /*
7514 * if the value is already no longer than the requested precision,
7515 * return the original value; similarly, if we don't have to do any
7516 * rounding to truncate to the requested precision, do not change
7517 * the original object; likewise, don't bother changing anything if
7518 * it's not a number
7519 */
7520 if (digits >= get_prec(ext) || get_dig(ext, digits) < 5
7521 || get_type(ext) != VMBN_T_NUM)
7522 {
7523 if (always_create)
7524 {
7525 /*
7526 * we must create a new object regardless, but it won't need
7527 * rounding
7528 */
7529 need_round = FALSE;
7530 }
7531 else
7532 {
7533 /* return the original value */
7534 new_val->set_nil();
7535 return ext;
7536 }
7537 }
7538
7539 /* allocate a new object with the requested precision */
7540 new_val->set_obj(create(vmg_ FALSE, digits));
7541 new_ext = get_objid_ext(vmg_ new_val->val.obj);
7542
7543 /* copy the sign, exponent, and type information */
7544 set_prec(new_ext, digits);
7545 set_neg(new_ext, get_neg(ext));
7546 set_exp(new_ext, get_exp(ext));
7547 set_type(new_ext, get_type(ext));
7548
7549 /*
7550 * if we don't need rounding, just truncate the old mantissa and
7551 * return the result
7552 */
7553 if (!need_round)
7554 {
7555 /* if the new size is smaller, truncate it */
7556 if (digits <= get_prec(ext))
7557 {
7558 /* copy the mantissa up to the requested new size */
7559 memcpy(new_ext + VMBN_MANT, ext + VMBN_MANT, (digits + 1)/2);
7560 }
7561 else
7562 {
7563 /* it's growing - simply copy the old value */
7564 copy_val(new_ext, ext, FALSE);
7565 }
7566
7567 /* return the new value */
7568 return new_ext;
7569 }
7570
7571 /* copy the mantissa up to the requested new precision */
7572 memcpy(new_ext + VMBN_MANT, ext + VMBN_MANT, (digits + 1)/2);
7573
7574 /* apply the rounding */
7575 for (need_carry = TRUE, idx = digits ; idx != 0 ; )
7576 {
7577 /* move to the previous index value */
7578 --idx;
7579
7580 /* round up - if it's a 9, we need to carry */
7581 if (get_dig(new_ext, idx) == 9)
7582 {
7583 /* make it a zero, and keep going to carry it */
7584 set_dig(new_ext, idx, 0);
7585 }
7586 else
7587 {
7588 /* bump it up */
7589 set_dig(new_ext, idx, get_dig(new_ext, idx) + 1);
7590
7591 /* no carrying required, so we can stop now */
7592 need_carry = FALSE;
7593 break;
7594 }
7595 }
7596
7597 /*
7598 * if we need to carry, we must have had all nines, in which case we
7599 * now have all zeroes - put a 1 in for the first digit, and
7600 * increase the exponent to account for the change
7601 */
7602 if (need_carry)
7603 {
7604 /* set the lead digit to 1 */
7605 set_dig(new_ext, 0, 1);
7606
7607 /* increase the exponent */
7608 set_exp(new_ext, get_exp(new_ext) + 1);
7609 }
7610
7611 /* return the new extension */
7612 return new_ext;
7613 }
7614
7615 /*
7616 * Convert a value to a big number value
7617 */
cvt_to_bignum(VMG_ vm_obj_id_t self,vm_val_t * val) const7618 int CVmObjBigNum::cvt_to_bignum(VMG_ vm_obj_id_t self, vm_val_t *val) const
7619 {
7620 /* if it's an integer, convert it to a BigNum value */
7621 if (val->typ == VM_INT)
7622 {
7623 /*
7624 * put my own value on the stack to ensure I'm not garbage
7625 * collected when creating the new object
7626 */
7627 G_stk->push()->set_obj(self);
7628
7629 /* it's an integer - convert it to a BigNum */
7630 val->set_obj(create(vmg_ FALSE, val->val.intval, 32));
7631
7632 /* done protecting my object reference */
7633 G_stk->discard();
7634 }
7635
7636 /* if it's not a BigNumberobject, we can't handle it */
7637 if (val->typ != VM_OBJ
7638 || (vm_objp(vmg_ val->val.obj)->get_metaclass_reg()
7639 != metaclass_reg_))
7640 {
7641 /* indicate that conversion was unsuccessful */
7642 return FALSE;
7643 }
7644
7645 /* successful conversion */
7646 return TRUE;
7647 }
7648
7649 /*
7650 * Allocate a temporary register
7651 */
alloc_temp_reg(VMG_ size_t prec,uint * hdl)7652 char *CVmObjBigNum::alloc_temp_reg(VMG_ size_t prec, uint *hdl)
7653 {
7654 char *p;
7655
7656 /* allocate a register with enough space for the desired precision */
7657 p = G_bignum_cache->alloc_reg(calc_alloc(prec), hdl);
7658
7659 /* if that succeeded, initialize the memory */
7660 if (p != 0)
7661 {
7662 /* set the desired precision */
7663 set_prec(p, prec);
7664
7665 /* initialize the flags */
7666 p[VMBN_FLAGS] = 0;
7667 }
7668
7669 /* return the register memory */
7670 return p;
7671 }
7672
7673 /*
7674 * release a temporary register
7675 */
release_temp_reg(VMG_ uint hdl)7676 void CVmObjBigNum::release_temp_reg(VMG_ uint hdl)
7677 {
7678 /* release the register to the cache */
7679 G_bignum_cache->release_reg(hdl);
7680 }
7681
7682 /*
7683 * Allocate a group of temporary registers
7684 */
alloc_temp_regs(VMG_ size_t prec,size_t cnt,...)7685 void CVmObjBigNum::alloc_temp_regs(VMG_ size_t prec, size_t cnt, ...)
7686 {
7687 va_list marker;
7688 size_t i;
7689 int failed;
7690 char **ext_ptr;
7691 uint *hdl_ptr;
7692
7693 /* set up to read varargs */
7694 va_start(marker, cnt);
7695
7696 /* no failures yet */
7697 failed = FALSE;
7698
7699 /* scan the varargs list */
7700 for (i = 0 ; i < cnt ; ++i)
7701 {
7702 /* get the next argument */
7703 ext_ptr = va_arg(marker, char **);
7704 hdl_ptr = va_arg(marker, uint *);
7705
7706 /* allocate a register */
7707 *ext_ptr = alloc_temp_reg(vmg_ prec, hdl_ptr);
7708
7709 /* if this allocation failed, note it, but keep going for now */
7710 if (*ext_ptr == 0)
7711 failed = TRUE;
7712 }
7713
7714 /* done reading argument */
7715 va_end(marker);
7716
7717 /* if we had any failures, free all of the registers we allocated */
7718 if (failed)
7719 {
7720 /* restart reading the varargs */
7721 va_start(marker, cnt);
7722
7723 /* scan the varargs and free the successfully allocated registers */
7724 for (i = 0 ; i < cnt ; ++i)
7725 {
7726 /* get the next argument */
7727 ext_ptr = va_arg(marker, char **);
7728 hdl_ptr = va_arg(marker, uint *);
7729
7730 /* free this register if we successfully allocated it */
7731 if (*ext_ptr != 0)
7732 release_temp_reg(vmg_ *hdl_ptr);
7733 }
7734
7735 /* done reading varargs */
7736 va_end(marker);
7737
7738 /* throw the error */
7739 err_throw(VMERR_BIGNUM_NO_REGS);
7740 }
7741 }
7742
7743 /*
7744 * Release a block of temporary registers
7745 */
release_temp_regs(VMG_ size_t cnt,...)7746 void CVmObjBigNum::release_temp_regs(VMG_ size_t cnt, ...)
7747 {
7748 size_t i;
7749 va_list marker;
7750
7751 /* start reading the varargs */
7752 va_start(marker, cnt);
7753
7754 /* scan the varargs and free the listed registers */
7755 for (i = 0 ; i < cnt ; ++i)
7756 {
7757 uint hdl;
7758
7759 /* get the next handle */
7760 hdl = va_arg(marker, uint);
7761
7762 /* free this register */
7763 release_temp_reg(vmg_ hdl);
7764 }
7765
7766 /* done reading varargs */
7767 va_end(marker);
7768 }
7769
7770
7771 /* ------------------------------------------------------------------------ */
7772 /*
7773 * Register cache
7774 */
7775
7776 /*
7777 * initialize
7778 */
CVmBigNumCache(size_t max_regs)7779 CVmBigNumCache::CVmBigNumCache(size_t max_regs)
7780 {
7781 CVmBigNumCacheReg *p;
7782 size_t i;
7783
7784 /* allocate our register array */
7785 reg_ = (CVmBigNumCacheReg *)t3malloc(max_regs * sizeof(reg_[0]));
7786
7787 /* remember the number of registers */
7788 max_regs_ = max_regs;
7789
7790 /* clear the list heads */
7791 free_reg_ = 0;
7792 unalloc_reg_ = 0;
7793
7794 /* we haven't actually allocated any registers yet - clear them out */
7795 for (p = reg_, i = max_regs ; i != 0 ; ++p, --i)
7796 {
7797 /* clear this register descriptor */
7798 p->clear();
7799
7800 /* link it into the unallocated list */
7801 p->nxt_ = unalloc_reg_;
7802 unalloc_reg_ = p;
7803 }
7804
7805 /* we haven't allocated the constants registers yet */
7806 ln10_.clear();
7807 pi_.clear();
7808 e_.clear();
7809 }
7810
7811 /*
7812 * delete
7813 */
~CVmBigNumCache()7814 CVmBigNumCache::~CVmBigNumCache()
7815 {
7816 CVmBigNumCacheReg *p;
7817 size_t i;
7818
7819 /* delete each of our allocated registers */
7820 for (p = reg_, i = max_regs_ ; i != 0 ; ++p, --i)
7821 p->free_mem();
7822
7823 /* free the register list array */
7824 t3free(reg_);
7825
7826 /* free the constant value registers */
7827 ln10_.free_mem();
7828 pi_.free_mem();
7829 e_.free_mem();
7830 }
7831
7832 /*
7833 * Allocate a register
7834 */
alloc_reg(size_t siz,uint * hdl)7835 char *CVmBigNumCache::alloc_reg(size_t siz, uint *hdl)
7836 {
7837 CVmBigNumCacheReg *p;
7838 CVmBigNumCacheReg *prv;
7839
7840 /*
7841 * search the free list for an available register satisfying the
7842 * size requirements
7843 */
7844 for (p = free_reg_, prv = 0 ; p != 0 ; prv = p, p = p->nxt_)
7845 {
7846 /* if it satisfies the size requirements, return it */
7847 if (p->siz_ >= siz)
7848 {
7849 /* unlink it from the free list */
7850 if (prv == 0)
7851 free_reg_ = p->nxt_;
7852 else
7853 prv->nxt_ = p->nxt_;
7854
7855 /* it's no longer in the free list */
7856 p->nxt_ = 0;
7857
7858 /* return it */
7859 *hdl = (uint)(p - reg_);
7860 return p->buf_;
7861 }
7862 }
7863
7864 /*
7865 * if there's an unallocated register, allocate it and use it;
7866 * otherwise, reallocate the smallest free register
7867 */
7868 if (unalloc_reg_ != 0)
7869 {
7870 /* use the first unallocated register */
7871 p = unalloc_reg_;
7872
7873 /* unlink it from the list */
7874 unalloc_reg_ = unalloc_reg_->nxt_;
7875 }
7876 else if (free_reg_ != 0)
7877 {
7878 CVmBigNumCacheReg *min_free_p;
7879 CVmBigNumCacheReg *min_free_prv;
7880
7881 /* search for the smallest free register */
7882 for (min_free_p = 0, p = free_reg_, prv = 0 ; p != 0 ;
7883 prv = p, p = p->nxt_)
7884 {
7885 /* if it's the smallest so far, remember it */
7886 if (min_free_p == 0 || p->siz_ < min_free_p->siz_)
7887 {
7888 /* remember it */
7889 min_free_p = p;
7890 min_free_prv = prv;
7891 }
7892 }
7893
7894 /* use the smallest register we found */
7895 p = min_free_p;
7896
7897 /* unlink it from the list */
7898 if (min_free_prv == 0)
7899 free_reg_ = p->nxt_;
7900 else
7901 min_free_prv->nxt_ = p->nxt_;
7902 }
7903 else
7904 {
7905 /* there are no free registers - return failure */
7906 return 0;
7907 }
7908
7909 /*
7910 * we found a register that was either previously unallocated, or
7911 * was previously allocated but was too small - allocate or
7912 * reallocate the register at the new size
7913 */
7914 p->alloc_mem(siz);
7915
7916 /* return the new register */
7917 *hdl = (uint)(p - reg_);
7918 return p->buf_;
7919 }
7920
7921 /*
7922 * Release a register
7923 */
release_reg(uint hdl)7924 void CVmBigNumCache::release_reg(uint hdl)
7925 {
7926 CVmBigNumCacheReg *p = ®_[hdl];
7927
7928 /* add the register to the free list */
7929 p->nxt_ = free_reg_;
7930 free_reg_ = p;
7931 }
7932
7933 /*
7934 * Release all registers
7935 */
release_all()7936 void CVmBigNumCache::release_all()
7937 {
7938 size_t i;
7939
7940 /* mark each of our registers as not in use */
7941 for (i = 0 ; i < max_regs_ ; ++i)
7942 release_reg(i);
7943 }
7944