1 /*
2 * Copyright (c) 2002 by The XFree86 Project, Inc.
3 *
4 * Permission is hereby granted, free of charge, to any person obtaining a
5 * copy of this software and associated documentation files (the "Software"),
6 * to deal in the Software without restriction, including without limitation
7 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
8 * and/or sell copies of the Software, and to permit persons to whom the
9 * Software is furnished to do so, subject to the following conditions:
10 *
11 * The above copyright notice and this permission notice shall be included in
12 * all copies or substantial portions of the Software.
13 *
14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
17 * THE XFREE86 PROJECT BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
18 * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF
19 * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
20 * SOFTWARE.
21 *
22 * Except as contained in this notice, the name of the XFree86 Project shall
23 * not be used in advertising or otherwise to promote the sale, use or other
24 * dealings in this Software without prior written authorization from the
25 * XFree86 Project.
26 *
27 * Author: Paulo César Pereira de Andrade
28 */
29
30 /* $XFree86: xc/programs/xedit/lisp/mathimp.c,v 1.14tsi Exp $ */
31
32
33 /*
34 * Defines
35 */
36 #ifdef __GNUC__
37 #define CONST __attribute__ ((__const__))
38 #else
39 #define CONST /**/
40 #endif
41
42 #define N_FIXNUM 1
43 #define N_BIGNUM 2
44 #define N_FLONUM 3
45 #define N_FIXRATIO 4
46 #define N_BIGRATIO 5
47
48 #define NOP_ADD 1
49 #define NOP_SUB 2
50 #define NOP_MUL 3
51 #define NOP_DIV 4
52
53 #define NDIVIDE_CEIL 1
54 #define NDIVIDE_FLOOR 2
55 #define NDIVIDE_ROUND 3
56 #define NDIVIDE_TRUNC 4
57
58 /* real part from number */
59 #define NREAL(num) &((num)->real)
60 #define NRTYPE(num) (num)->real.type
61 #define NRFI(num) (num)->real.data.fixnum
62 #define NRBI(num) (num)->real.data.bignum
63 #define NRFF(num) (num)->real.data.flonum
64 #define NRFRN(Num) (Num)->real.data.fixratio.num
65 #define NRFRD(num) (num)->real.data.fixratio.den
66 #define NRBR(num) (num)->real.data.bigratio
67 #define NRBRN(num) mpr_num(NRBR(num))
68 #define NRBRD(num) mpr_den(NRBR(num))
69
70 #define NRCLEAR_BI(num) mpi_clear(NRBI(num)); XFREE(NRBI(num))
71 #define NRCLEAR_BR(num) mpr_clear(NRBR(num)); XFREE(NRBR(num))
72
73 /* imag part from number */
74 #define NIMAG(num) &((num)->imag)
75 #define NITYPE(num) (num)->imag.type
76 #define NIFI(num) (num)->imag.data.fixnum
77 #define NIBI(num) (num)->imag.data.bignum
78 #define NIFF(num) (num)->imag.data.flonum
79 #define NIFRN(Num) (Num)->imag.data.fixratio.num
80 #define NIFRD(num) (num)->imag.data.fixratio.den
81 #define NIBR(num) (num)->imag.data.bigratio
82 #define NIBRN(obj) mpr_num(NIBR(obj))
83 #define NIBRD(obj) mpr_den(NIBR(obj))
84
85 /* real number fields */
86 #define RTYPE(real) (real)->type
87 #define RFI(real) (real)->data.fixnum
88 #define RBI(real) (real)->data.bignum
89 #define RFF(real) (real)->data.flonum
90 #define RFRN(real) (real)->data.fixratio.num
91 #define RFRD(real) (real)->data.fixratio.den
92 #define RBR(real) (real)->data.bigratio
93 #define RBRN(real) mpr_num(RBR(real))
94 #define RBRD(real) mpr_den(RBR(real))
95
96 #define RINTEGERP(real) \
97 (RTYPE(real) == N_FIXNUM || RTYPE(real) == N_BIGNUM)
98
99 #define RCLEAR_BI(real) mpi_clear(RBI(real)); XFREE(RBI(real))
100 #define RCLEAR_BR(real) mpr_clear(RBR(real)); XFREE(RBR(real))
101
102 /* numeric value from lisp object */
103 #define OFI(object) FIXNUM_VALUE(object)
104 #define OII(object) INT_VALUE(object)
105 #define OBI(object) (object)->data.mp.integer
106 #define ODF(object) DFLOAT_VALUE(object)
107 #define OFRN(object) (object)->data.ratio.numerator
108 #define OFRD(object) (object)->data.ratio.denominator
109 #define OBR(object) (object)->data.mp.ratio
110 #define OBRN(object) mpr_num(OBR(object))
111 #define OBRD(object) mpr_den(OBR(object))
112 #define OCXR(object) (object)->data.complex.real
113 #define OCXI(object) (object)->data.complex.imag
114
115 #define XALLOC(type) LispMalloc(sizeof(type))
116 #define XFREE(ptr) LispFree(ptr)
117
118
119 /*
120 * Types
121 */
122 typedef struct _n_real {
123 char type;
124 union {
125 long fixnum;
126 mpi *bignum;
127 double flonum;
128 struct {
129 long num;
130 long den;
131 } fixratio;
132 mpr *bigratio;
133 } data;
134 } n_real;
135
136 typedef struct _n_number {
137 char complex;
138 n_real real;
139 n_real imag;
140 } n_number;
141
142
143 /*
144 * Prototypes
145 */
146 static void number_init(void);
147 static LispObj *number_pi(void);
148
149 static void set_real_real(n_real*, n_real*);
150 static void set_real_object(n_real*, LispObj*);
151 static void set_number_object(n_number*, LispObj*);
152 static void clear_real(n_real*);
153 static void clear_number(n_number*);
154
155 static LispObj *make_real_object(n_real*);
156 static LispObj *make_number_object(n_number*);
157
158 static void fatal_error(int);
159 static void fatal_object_error(LispObj*, int);
160 static void fatal_builtin_object_error(LispBuiltin*, LispObj*, int);
161
162 static double bi_getd(mpi*);
163 static double br_getd(mpr*);
164
165 /* add */
166 static void add_real_object(n_real*, LispObj*);
167 static void add_number_object(n_number*, LispObj*);
168
169 /* sub */
170 static void sub_real_object(n_real*, LispObj*);
171 static void sub_number_object(n_number*, LispObj*);
172
173 /* mul */
174 static void mul_real_object(n_real*, LispObj*);
175 static void mul_number_object(n_number*, LispObj*);
176
177 /* div */
178 static void div_real_object(n_real*, LispObj*);
179 static void div_number_object(n_number*, LispObj*);
180
181 /* compare */
182 static int cmp_real_real(n_real*, n_real*);
183 static int cmp_real_object(n_real*, LispObj*);
184 #if 0 /* not used */
185 static int cmp_number_object(n_number*, LispObj*);
186 #endif
187 static int cmp_object_object(LispObj*, LispObj*, int);
188
189 /* fixnum */
190 static INLINE int fi_fi_add_overflow(long, long) CONST;
191 static INLINE int fi_fi_sub_overflow(long, long) CONST;
192 static INLINE int fi_fi_mul_overflow(long, long) CONST;
193
194 /* bignum */
195 static void rbi_canonicalize(n_real*);
196
197 /* ratio */
198 static void rfr_canonicalize(n_real*);
199 static void rbr_canonicalize(n_real*);
200
201 /* complex */
202 static void ncx_canonicalize(n_number*);
203
204 /* abs */
205 static void abs_real(n_real*);
206 static void abs_number(n_number*);
207 static void nabs_cx(n_number*);
208 static INLINE void rabs_fi(n_real*);
209 static INLINE void rabs_bi(n_real*);
210 static INLINE void rabs_ff(n_real*);
211 static INLINE void rabs_fr(n_real*);
212 static INLINE void rabs_br(n_real*);
213
214 /* neg */
215 static void neg_real(n_real*);
216 static void neg_number(n_number*);
217 static void rneg_fi(n_real*);
218 static INLINE void rneg_bi(n_real*);
219 static INLINE void rneg_ff(n_real*);
220 static INLINE void rneg_fr(n_real*);
221 static INLINE void rneg_br(n_real*);
222
223 /* sqrt */
224 static void sqrt_real(n_real*);
225 static void sqrt_number(n_number*);
226 static void rsqrt_xi(n_real*);
227 static void rsqrt_xr(n_real*);
228 static void rsqrt_ff(n_real*);
229 static void nsqrt_cx(n_number*);
230 static void nsqrt_xi(n_number*);
231 static void nsqrt_ff(n_number*);
232 static void nsqrt_xr(n_number*);
233
234 /* mod */
235 static void mod_real_real(n_real*, n_real*);
236 static void mod_real_object(n_real*, LispObj*);
237 static void rmod_fi_fi(n_real*, long);
238 static void rmod_fi_bi(n_real*, mpi*);
239 static void rmod_bi_fi(n_real*, long);
240 static void rmod_bi_bi(n_real*, mpi*);
241
242 /* rem */
243 static void rem_real_object(n_real*, LispObj*);
244 static void rrem_fi_fi(n_real*, long);
245 static void rrem_fi_bi(n_real*, mpi*);
246 static void rrem_bi_fi(n_real*, long);
247 static void rrem_bi_bi(n_real*, mpi*);
248
249 /* gcd */
250 static void gcd_real_object(n_real*, LispObj*);
251
252 /* and */
253 static void and_real_object(n_real*, LispObj*);
254
255 /* eqv */
256 static void eqv_real_object(n_real*, LispObj*);
257
258 /* ior */
259 static void ior_real_object(n_real*, LispObj*);
260
261 /* not */
262 static void not_real(n_real*);
263
264 /* xor */
265 static void xor_real_object(n_real*, LispObj*);
266
267 /* divide */
268 static void divide_number_object(n_number*, LispObj*, int, int);
269 static void ndivide_xi_xi(n_number*, LispObj*, int, int);
270 static void ndivide_flonum(n_number*, double, double, int, int);
271 static void ndivide_xi_xr(n_number*, LispObj*, int, int);
272 static void ndivide_xr_xi(n_number*, LispObj*, int, int);
273 static void ndivide_xr_xr(n_number*, LispObj*, int, int);
274
275 /* real complex */
276 static void nadd_re_cx(n_number*, LispObj*);
277 static void nsub_re_cx(n_number*, LispObj*);
278 static void nmul_re_cx(n_number*, LispObj*);
279 static void ndiv_re_cx(n_number*, LispObj*);
280
281 /* complex real */
282 static void nadd_cx_re(n_number*, LispObj*);
283 static void nsub_cx_re(n_number*, LispObj*);
284 static void nmul_cx_re(n_number*, LispObj*);
285 static void ndiv_cx_re(n_number*, LispObj*);
286
287 /* complex complex */
288 static void nadd_cx_cx(n_number*, LispObj*);
289 static void nsub_cx_cx(n_number*, LispObj*);
290 static void nmul_cx_cx(n_number*, LispObj*);
291 static void ndiv_cx_cx(n_number*, LispObj*);
292 static int cmp_cx_cx(LispObj*, LispObj*);
293
294 /* flonum flonum */
295 static void radd_flonum(n_real*, double, double);
296 static void rsub_flonum(n_real*, double, double);
297 static void rmul_flonum(n_real*, double, double);
298 static void rdiv_flonum(n_real*, double, double);
299 static int cmp_flonum(double, double);
300
301 /* fixnum fixnum */
302 static void rop_fi_fi_bi(n_real*, long, int);
303 static INLINE void radd_fi_fi(n_real*, long);
304 static INLINE void rsub_fi_fi(n_real*, long);
305 static INLINE void rmul_fi_fi(n_real*, long);
306 static INLINE void rdiv_fi_fi(n_real*, long);
307 static INLINE int cmp_fi_fi(long, long);
308 static void ndivide_fi_fi(n_number*, long, int, int);
309
310 /* fixnum bignum */
311 static void rop_fi_bi_xi(n_real*, mpi*, int);
312 static INLINE void radd_fi_bi(n_real*, mpi*);
313 static INLINE void rsub_fi_bi(n_real*, mpi*);
314 static INLINE void rmul_fi_bi(n_real*, mpi*);
315 static void rdiv_fi_bi(n_real*, mpi*);
316 static INLINE int cmp_fi_bi(long, mpi*);
317
318 /* fixnum fixratio */
319 static void rop_fi_fr_as_xr(n_real*, long, long, int);
320 static void rop_fi_fr_md_xr(n_real*, long, long, int);
321 static INLINE void radd_fi_fr(n_real*, long, long);
322 static INLINE void rsub_fi_fr(n_real*, long, long);
323 static INLINE void rmul_fi_fr(n_real*, long, long);
324 static INLINE void rdiv_fi_fr(n_real*, long, long);
325 static INLINE int cmp_fi_fr(long, long, long);
326
327 /* fixnum bigratio */
328 static void rop_fi_br_as_xr(n_real*, mpr*, int);
329 static void rop_fi_br_md_xr(n_real*, mpr*, int);
330 static INLINE void radd_fi_br(n_real*, mpr*);
331 static INLINE void rsub_fi_br(n_real*, mpr*);
332 static INLINE void rmul_fi_br(n_real*, mpr*);
333 static INLINE void rdiv_fi_br(n_real*, mpr*);
334 static INLINE int cmp_fi_br(long, mpr*);
335
336 /* bignum fixnum */
337 static INLINE void radd_bi_fi(n_real*, long);
338 static INLINE void rsub_bi_fi(n_real*, long);
339 static INLINE void rmul_bi_fi(n_real*, long);
340 static void rdiv_bi_fi(n_real*, long);
341 static INLINE int cmp_bi_fi(mpi*, long);
342
343 /* bignum bignum */
344 static INLINE void radd_bi_bi(n_real*, mpi*);
345 static INLINE void rsub_bi_bi(n_real*, mpi*);
346 static INLINE void rmul_bi_bi(n_real*, mpi*);
347 static void rdiv_bi_bi(n_real*, mpi*);
348 static INLINE int cmp_bi_bi(mpi*, mpi*);
349
350 /* bignum fixratio */
351 static void rop_bi_fr_as_xr(n_real*, long, long, int);
352 static void rop_bi_fr_md_xr(n_real*, long, long, int);
353 static INLINE void radd_bi_fr(n_real*, long, long);
354 static INLINE void rsub_bi_fr(n_real*, long, long);
355 static INLINE void rmul_bi_fr(n_real*, long, long);
356 static INLINE void rdiv_bi_fr(n_real*, long, long);
357 static int cmp_bi_fr(mpi*, long, long);
358
359 /* bignum bigratio */
360 static void rop_bi_br_as_xr(n_real*, mpr*, int);
361 static void rop_bi_br_md_xr(n_real*, mpr*, int);
362 static INLINE void radd_bi_br(n_real*, mpr*);
363 static INLINE void rsub_bi_br(n_real*, mpr*);
364 static INLINE void rmul_bi_br(n_real*, mpr*);
365 static INLINE void rdiv_bi_br(n_real*, mpr*);
366 static int cmp_bi_br(mpi*, mpr*);
367
368 /* fixratio fixnum */
369 static void rop_fr_fi_as_xr(n_real*, long, int);
370 static void rop_fr_fi_md_xr(n_real*, long, int);
371 static INLINE void radd_fr_fi(n_real*, long);
372 static INLINE void rsub_fr_fi(n_real*, long);
373 static INLINE void rmul_fr_fi(n_real*, long);
374 static INLINE void rdiv_fr_fi(n_real*, long);
375 static INLINE int cmp_fr_fi(long, long, long);
376
377 /* fixratio bignum */
378 static void rop_fr_bi_as_xr(n_real*, mpi*, int);
379 static void rop_fr_bi_md_xr(n_real*, mpi*, int);
380 static INLINE void radd_fr_bi(n_real*, mpi*);
381 static INLINE void rsub_fr_bi(n_real*, mpi*);
382 static INLINE void rmul_fr_bi(n_real*, mpi*);
383 static INLINE void rdiv_fr_bi(n_real*, mpi*);
384 static int cmp_fr_bi(long, long, mpi*);
385
386 /* fixratio fixratio */
387 static void rop_fr_fr_as_xr(n_real*, long, long, int);
388 static void rop_fr_fr_md_xr(n_real*, long, long, int);
389 static INLINE void radd_fr_fr(n_real*, long, long);
390 static INLINE void rsub_fr_fr(n_real*, long, long);
391 static INLINE void rmul_fr_fr(n_real*, long, long);
392 static INLINE void rdiv_fr_fr(n_real*, long, long);
393 static INLINE int cmp_fr_fr(long, long, long, long);
394
395 /* fixratio bigratio */
396 static void rop_fr_br_asmd_xr(n_real*, mpr*, int);
397 static INLINE void radd_fr_br(n_real*, mpr*);
398 static INLINE void rsub_fr_br(n_real*, mpr*);
399 static INLINE void rmul_fr_br(n_real*, mpr*);
400 static INLINE void rdiv_fr_br(n_real*, mpr*);
401 static int cmp_fr_br(long, long, mpr*);
402
403 /* bigratio fixnum */
404 static void rop_br_fi_asmd_xr(n_real*, long, int);
405 static INLINE void radd_br_fi(n_real*, long);
406 static INLINE void rsub_br_fi(n_real*, long);
407 static INLINE void rmul_br_fi(n_real*, long);
408 static INLINE void rdiv_br_fi(n_real*, long);
409 static int cmp_br_fi(mpr*, long);
410
411 /* bigratio bignum */
412 static void rop_br_bi_as_xr(n_real*, mpi*, int);
413 static INLINE void radd_br_bi(n_real*, mpi*);
414 static INLINE void rsub_br_bi(n_real*, mpi*);
415 static INLINE void rmul_br_bi(n_real*, mpi*);
416 static INLINE void rdiv_br_bi(n_real*, mpi*);
417 static int cmp_br_bi(mpr*, mpi*);
418
419 /* bigratio fixratio */
420 static void rop_br_fr_asmd_xr(n_real*, long, long, int);
421 static INLINE void radd_br_fr(n_real*, long, long);
422 static INLINE void rsub_br_fr(n_real*, long, long);
423 static INLINE void rmul_br_fr(n_real*, long, long);
424 static INLINE void rdiv_br_fr(n_real*, long, long);
425 static int cmp_br_fr(mpr*, long, long);
426
427 /* bigratio bigratio */
428 static INLINE void radd_br_br(n_real*, mpr*);
429 static INLINE void rsub_br_br(n_real*, mpr*);
430 static INLINE void rmul_br_br(n_real*, mpr*);
431 static INLINE void rdiv_br_br(n_real*, mpr*);
432 static INLINE int cmp_br_br(mpr*, mpr*);
433
434 /*
435 * Initialization
436 */
437 static n_real zero, one, two;
438
439 static const char *fatal_error_strings[] = {
440 #define DIVIDE_BY_ZERO 0
441 "divide by zero",
442 #define FLOATING_POINT_OVERFLOW 1
443 "floating point overflow",
444 #define FLOATING_POINT_EXCEPTION 2
445 "floating point exception"
446 };
447
448 static const char *fatal_object_error_strings[] = {
449 #define NOT_A_NUMBER 0
450 "is not a number",
451 #define NOT_A_REAL_NUMBER 1
452 "is not a real number",
453 #define NOT_AN_INTEGER 2
454 "is not an integer"
455 };
456
457 /*
458 * Implementation
459 */
460 static void
fatal_error(int num)461 fatal_error(int num)
462 {
463 LispDestroy("%s", fatal_error_strings[num]);
464 }
465
466 static void
fatal_object_error(LispObj * obj,int num)467 fatal_object_error(LispObj *obj, int num)
468 {
469 LispDestroy("%s %s", STROBJ(obj), fatal_object_error_strings[num]);
470 }
471
472 static void
fatal_builtin_object_error(LispBuiltin * builtin,LispObj * obj,int num)473 fatal_builtin_object_error(LispBuiltin *builtin, LispObj *obj, int num)
474 {
475 LispDestroy("%s: %s %s", STRFUN(builtin), STROBJ(obj),
476 fatal_object_error_strings[num]);
477 }
478
479 static void
number_init(void)480 number_init(void)
481 {
482 zero.type = one.type = two.type = N_FIXNUM;
483 zero.data.fixnum = 0;
484 one.data.fixnum = 1;
485 two.data.fixnum = 2;
486 }
487
488 static double
bi_getd(mpi * bignum)489 bi_getd(mpi *bignum)
490 {
491 double value = mpi_getd(bignum);
492
493 if (!finite(value))
494 fatal_error(FLOATING_POINT_EXCEPTION);
495
496 return (value);
497 }
498
499 static double
br_getd(mpr * bigratio)500 br_getd(mpr *bigratio)
501 {
502 double value = mpr_getd(bigratio);
503
504 if (!finite(value))
505 fatal_error(FLOATING_POINT_EXCEPTION);
506
507 return (value);
508 }
509
510 static LispObj *
number_pi(void)511 number_pi(void)
512 {
513 LispObj *result;
514 #ifndef M_PI
515 #define M_PI 3.14159265358979323846
516 #endif
517 result = DFLOAT(M_PI);
518
519 return (result);
520 }
521
522 static void
set_real_real(n_real * real,n_real * val)523 set_real_real(n_real *real, n_real *val)
524 {
525 switch (RTYPE(real) = RTYPE(val)) {
526 case N_FIXNUM:
527 RFI(real) = RFI(val);
528 break;
529 case N_BIGNUM:
530 RBI(real) = XALLOC(mpi);
531 mpi_init(RBI(real));
532 mpi_set(RBI(real), RBI(val));
533 break;
534 case N_FLONUM:
535 RFF(real) = RFF(val);
536 break;
537 case N_FIXRATIO:
538 RFRN(real) = RFRN(val);
539 RFRD(real) = RFRD(val);
540 break;
541 case N_BIGRATIO:
542 RBR(real) = XALLOC(mpr);
543 mpr_init(RBR(real));
544 mpr_set(RBR(real), RBR(val));
545 break;
546 }
547 }
548
549 static void
set_real_object(n_real * real,LispObj * obj)550 set_real_object(n_real *real, LispObj *obj)
551 {
552 switch (OBJECT_TYPE(obj)) {
553 case LispFixnum_t:
554 RTYPE(real) = N_FIXNUM;
555 RFI(real) = OFI(obj);
556 break;
557 case LispInteger_t:
558 RTYPE(real) = N_FIXNUM;
559 RFI(real) = OII(obj);
560 break;
561 case LispBignum_t:
562 RTYPE(real) = N_BIGNUM;
563 RBI(real) = XALLOC(mpi);
564 mpi_init(RBI(real));
565 mpi_set(RBI(real), OBI(obj));
566 break;
567 case LispDFloat_t:
568 RTYPE(real) = N_FLONUM;
569 RFF(real) = ODF(obj);
570 break;
571 case LispRatio_t:
572 RTYPE(real) = N_FIXRATIO;
573 RFRN(real) = OFRN(obj);
574 RFRD(real) = OFRD(obj);
575 break;
576 case LispBigratio_t:
577 RTYPE(real) = N_BIGRATIO;
578 RBR(real) = XALLOC(mpr);
579 mpr_init(RBR(real));
580 mpr_set(RBR(real), OBR(obj));
581 break;
582 default:
583 fatal_object_error(obj, NOT_A_REAL_NUMBER);
584 break;
585 }
586 }
587
588 static void
set_number_object(n_number * num,LispObj * obj)589 set_number_object(n_number *num, LispObj *obj)
590 {
591 switch (OBJECT_TYPE(obj)) {
592 case LispFixnum_t:
593 num->complex = 0;
594 NRTYPE(num) = N_FIXNUM;
595 NRFI(num) = OFI(obj);
596 break;
597 case LispInteger_t:
598 num->complex = 0;
599 NRTYPE(num) = N_FIXNUM;
600 NRFI(num) = OII(obj);
601 break;
602 case LispBignum_t:
603 num->complex = 0;
604 NRTYPE(num) = N_BIGNUM;
605 NRBI(num) = XALLOC(mpi);
606 mpi_init(NRBI(num));
607 mpi_set(NRBI(num), OBI(obj));
608 break;
609 case LispDFloat_t:
610 num->complex = 0;
611 NRTYPE(num) = N_FLONUM;
612 NRFF(num) = ODF(obj);
613 break;
614 case LispRatio_t:
615 num->complex = 0;
616 NRTYPE(num) = N_FIXRATIO;
617 NRFRN(num) = OFRN(obj);
618 NRFRD(num) = OFRD(obj);
619 break;
620 case LispBigratio_t:
621 num->complex = 0;
622 NRTYPE(num) = N_BIGRATIO;
623 NRBR(num) = XALLOC(mpr);
624 mpr_init(NRBR(num));
625 mpr_set(NRBR(num), OBR(obj));
626 break;
627 case LispComplex_t:
628 num->complex = 1;
629 set_real_object(NREAL(num), OCXR(obj));
630 set_real_object(NIMAG(num), OCXI(obj));
631 break;
632 default:
633 fatal_object_error(obj, NOT_A_NUMBER);
634 break;
635 }
636 }
637
638 static void
clear_real(n_real * real)639 clear_real(n_real *real)
640 {
641 if (RTYPE(real) == N_BIGNUM) {
642 mpi_clear(RBI(real));
643 XFREE(RBI(real));
644 }
645 else if (RTYPE(real) == N_BIGRATIO) {
646 mpr_clear(RBR(real));
647 XFREE(RBR(real));
648 }
649 }
650
651 static void
clear_number(n_number * num)652 clear_number(n_number *num)
653 {
654 clear_real(NREAL(num));
655 if (num->complex)
656 clear_real(NIMAG(num));
657 }
658
659 static LispObj *
make_real_object(n_real * real)660 make_real_object(n_real *real)
661 {
662 LispObj *obj;
663
664 switch (RTYPE(real)) {
665 case N_FIXNUM:
666 if (RFI(real) > MOST_POSITIVE_FIXNUM ||
667 RFI(real) < MOST_NEGATIVE_FIXNUM) {
668 obj = LispNew(NIL, NIL);
669 obj->type = LispInteger_t;
670 OII(obj) = RFI(real);
671 }
672 else
673 obj = FIXNUM(RFI(real));
674 break;
675 case N_BIGNUM:
676 obj = BIGNUM(RBI(real));
677 break;
678 case N_FLONUM:
679 obj = DFLOAT(RFF(real));
680 break;
681 case N_FIXRATIO:
682 obj = LispNew(NIL, NIL);
683 obj->type = LispRatio_t;
684 OFRN(obj) = RFRN(real);
685 OFRD(obj) = RFRD(real);
686 break;
687 case N_BIGRATIO:
688 obj = BIGRATIO(RBR(real));
689 break;
690 default:
691 obj = NIL;
692 break;
693 }
694
695 return (obj);
696 }
697
698 static LispObj *
make_number_object(n_number * num)699 make_number_object(n_number *num)
700 {
701 LispObj *obj;
702
703 if (num->complex) {
704 GC_ENTER();
705
706 obj = LispNew(NIL, NIL);
707 GC_PROTECT(obj);
708 OCXI(obj) = NIL;
709 obj->type = LispComplex_t;
710 OCXR(obj) = make_real_object(NREAL(num));
711 OCXI(obj) = make_real_object(NIMAG(num));
712 GC_LEAVE();
713 }
714 else {
715 switch (NRTYPE(num)) {
716 case N_FIXNUM:
717 if (NRFI(num) > MOST_POSITIVE_FIXNUM ||
718 NRFI(num) < MOST_NEGATIVE_FIXNUM) {
719 obj = LispNew(NIL, NIL);
720 obj->type = LispInteger_t;
721 OII(obj) = NRFI(num);
722 }
723 else
724 obj = FIXNUM(NRFI(num));
725 break;
726 case N_BIGNUM:
727 obj = BIGNUM(NRBI(num));
728 break;
729 case N_FLONUM:
730 obj = DFLOAT(NRFF(num));
731 break;
732 case N_FIXRATIO:
733 obj = LispNew(NIL, NIL);
734 obj->type = LispRatio_t;
735 OFRN(obj) = NRFRN(num);
736 OFRD(obj) = NRFRD(num);
737 break;
738 case N_BIGRATIO:
739 obj = BIGRATIO(NRBR(num));
740 break;
741 default:
742 obj = NIL;
743 break;
744 }
745 }
746
747 return (obj);
748 }
749
750 #define DEFOP_REAL_REAL(OP) \
751 OP##_real_real(n_real *real, n_real *val) \
752 { \
753 switch (RTYPE(real)) { \
754 case N_FIXNUM: \
755 switch (RTYPE(val)) { \
756 case N_FIXNUM: \
757 r##OP##_fi_fi(real, RFI(val)); \
758 break; \
759 case N_BIGNUM: \
760 r##OP##_fi_bi(real, RBI(val)); \
761 break; \
762 case N_FLONUM: \
763 r##OP##_flonum(real, (double)RFI(real), RFF(val)); \
764 break; \
765 case N_FIXRATIO: \
766 r##OP##_fi_fr(real, RFRN(val), RFRD(val)); \
767 break; \
768 case N_BIGRATIO: \
769 r##OP##_fi_br(real, RBR(val)); \
770 break; \
771 } \
772 break; \
773 case N_BIGNUM: \
774 switch (RTYPE(val)) { \
775 case N_FIXNUM: \
776 r##OP##_bi_fi(real, RFI(val)); \
777 break; \
778 case N_BIGNUM: \
779 r##OP##_bi_bi(real, RBI(val)); \
780 break; \
781 case N_FLONUM: \
782 r##OP##_flonum(real, bi_getd(RBI(real)), RFF(val)); \
783 break; \
784 case N_FIXRATIO: \
785 r##OP##_bi_fr(real, RFRN(val), RFRD(val)); \
786 break; \
787 case N_BIGRATIO: \
788 r##OP##_bi_br(real, RBR(val)); \
789 break; \
790 } \
791 break; \
792 case N_FLONUM: \
793 switch (RTYPE(val)) { \
794 case N_FIXNUM: \
795 r##OP##_flonum(real, RFF(real), (double)RFI(val)); \
796 break; \
797 case N_BIGNUM: \
798 r##OP##_flonum(real, RFF(real), bi_getd(RBI(val))); \
799 break; \
800 case N_FLONUM: \
801 r##OP##_flonum(real, RFF(real), RFF(val)); \
802 break; \
803 case N_FIXRATIO: \
804 r##OP##_flonum(real, RFF(real), \
805 (double)RFRN(val) / (double)RFRD(val));\
806 break; \
807 case N_BIGRATIO: \
808 r##OP##_flonum(real, RFF(real), br_getd(RBR(val))); \
809 break; \
810 } \
811 break; \
812 case N_FIXRATIO: \
813 switch (RTYPE(val)) { \
814 case N_FIXNUM: \
815 r##OP##_fr_fi(real, RFI(val)); \
816 break; \
817 case N_BIGNUM: \
818 r##OP##_fr_bi(real, RBI(val)); \
819 break; \
820 case N_FLONUM: \
821 r##OP##_flonum(real, \
822 (double)RFRN(real) / (double)RFRD(real),\
823 RFF(val)); \
824 break; \
825 case N_FIXRATIO: \
826 r##OP##_fr_fr(real, RFRN(val), RFRD(val)); \
827 break; \
828 case N_BIGRATIO: \
829 r##OP##_fr_br(real, RBR(val)); \
830 break; \
831 } \
832 break; \
833 case N_BIGRATIO: \
834 switch (RTYPE(val)) { \
835 case N_FIXNUM: \
836 r##OP##_br_fi(real, RFI(val)); \
837 break; \
838 case N_BIGNUM: \
839 r##OP##_br_bi(real, RBI(val)); \
840 break; \
841 case N_FLONUM: \
842 r##OP##_flonum(real, br_getd(RBR(real)), RFF(val)); \
843 break; \
844 case N_FIXRATIO: \
845 r##OP##_br_fr(real, RFRN(val), RFRD(val)); \
846 break; \
847 case N_BIGRATIO: \
848 r##OP##_br_br(real, RBR(val)); \
849 break; \
850 } \
851 break; \
852 } \
853 }
854
855 static void
DEFOP_REAL_REAL(add)856 DEFOP_REAL_REAL(add)
857
858 static void
859 DEFOP_REAL_REAL(sub)
860
861 static void
862 DEFOP_REAL_REAL(div)
863
864 static void
865 DEFOP_REAL_REAL(mul)
866
867
868 #define DEFOP_REAL_OBJECT(OP) \
869 OP##_real_object(n_real *real, LispObj *obj) \
870 { \
871 switch (OBJECT_TYPE(obj)) { \
872 case LispFixnum_t: \
873 switch (RTYPE(real)) { \
874 case N_FIXNUM: \
875 r##OP##_fi_fi(real, OFI(obj)); \
876 break; \
877 case N_BIGNUM: \
878 r##OP##_bi_fi(real, OFI(obj)); \
879 break; \
880 case N_FLONUM: \
881 r##OP##_flonum(real, RFF(real), (double)OFI(obj)); \
882 break; \
883 case N_FIXRATIO: \
884 r##OP##_fr_fi(real, OFI(obj)); \
885 break; \
886 case N_BIGRATIO: \
887 r##OP##_br_fi(real, OFI(obj)); \
888 break; \
889 } \
890 break; \
891 case LispInteger_t: \
892 switch (RTYPE(real)) { \
893 case N_FIXNUM: \
894 r##OP##_fi_fi(real, OII(obj)); \
895 break; \
896 case N_BIGNUM: \
897 r##OP##_bi_fi(real, OII(obj)); \
898 break; \
899 case N_FLONUM: \
900 r##OP##_flonum(real, RFF(real), (double)OII(obj)); \
901 break; \
902 case N_FIXRATIO: \
903 r##OP##_fr_fi(real, OII(obj)); \
904 break; \
905 case N_BIGRATIO: \
906 r##OP##_br_fi(real, OII(obj)); \
907 break; \
908 } \
909 break; \
910 case LispBignum_t: \
911 switch (RTYPE(real)) { \
912 case N_FIXNUM: \
913 r##OP##_fi_bi(real, OBI(obj)); \
914 break; \
915 case N_BIGNUM: \
916 r##OP##_bi_bi(real, OBI(obj)); \
917 break; \
918 case N_FLONUM: \
919 r##OP##_flonum(real, RFF(real), bi_getd(OBI(obj))); \
920 break; \
921 case N_FIXRATIO: \
922 r##OP##_fr_bi(real, OBI(obj)); \
923 break; \
924 case N_BIGRATIO: \
925 r##OP##_br_bi(real, OBI(obj)); \
926 break; \
927 } \
928 break; \
929 case LispDFloat_t: \
930 switch (RTYPE(real)) { \
931 case N_FIXNUM: \
932 r##OP##_flonum(real, (double)RFI(real), ODF(obj)); \
933 break; \
934 case N_BIGNUM: \
935 r##OP##_flonum(real, bi_getd(RBI(real)), ODF(obj)); \
936 break; \
937 case N_FLONUM: \
938 r##OP##_flonum(real, RFF(real), ODF(obj)); \
939 break; \
940 case N_FIXRATIO: \
941 r##OP##_flonum(real, \
942 (double)RFRN(real) / (double)RFRD(real),\
943 ODF(obj)); \
944 break; \
945 case N_BIGRATIO: \
946 r##OP##_flonum(real, br_getd(RBR(real)), ODF(obj)); \
947 break; \
948 } \
949 break; \
950 case LispRatio_t: \
951 switch (RTYPE(real)) { \
952 case N_FIXNUM: \
953 r##OP##_fi_fr(real, OFRN(obj), OFRD(obj)); \
954 break; \
955 case N_BIGNUM: \
956 r##OP##_bi_fr(real, OFRN(obj), OFRD(obj)); \
957 break; \
958 case N_FLONUM: \
959 r##OP##_flonum(real, RFF(real), \
960 (double)OFRN(obj) / (double)OFRD(obj)); \
961 break; \
962 case N_FIXRATIO: \
963 r##OP##_fr_fr(real, OFRN(obj), OFRD(obj)); \
964 break; \
965 case N_BIGRATIO: \
966 r##OP##_br_fr(real, OFRN(obj), OFRD(obj)); \
967 break; \
968 } \
969 break; \
970 case LispBigratio_t: \
971 switch (RTYPE(real)) { \
972 case N_FIXNUM: \
973 r##OP##_fi_br(real, OBR(obj)); \
974 break; \
975 case N_BIGNUM: \
976 r##OP##_bi_br(real, OBR(obj)); \
977 break; \
978 case N_FLONUM: \
979 r##OP##_flonum(real, RFF(real), br_getd(OBR(obj))); \
980 break; \
981 case N_FIXRATIO: \
982 r##OP##_fr_br(real, OBR(obj)); \
983 break; \
984 case N_BIGRATIO: \
985 r##OP##_br_br(real, OBR(obj)); \
986 break; \
987 } \
988 break; \
989 default: \
990 fatal_object_error(obj, NOT_A_REAL_NUMBER); \
991 break; \
992 } \
993 }
994
995 static void
996 DEFOP_REAL_OBJECT(add)
997
998 static void
999 DEFOP_REAL_OBJECT(sub)
1000
1001 static void
1002 DEFOP_REAL_OBJECT(div)
1003
1004 static void
1005 DEFOP_REAL_OBJECT(mul)
1006
1007
1008 #define DEFOP_NUMBER_OBJECT(OP) \
1009 OP##_number_object(n_number *num, LispObj *obj) \
1010 { \
1011 if (num->complex) { \
1012 switch (OBJECT_TYPE(obj)) { \
1013 case LispFixnum_t: \
1014 case LispInteger_t: \
1015 case LispBignum_t: \
1016 case LispDFloat_t: \
1017 case LispRatio_t: \
1018 case LispBigratio_t: \
1019 n##OP##_cx_re(num, obj); \
1020 break; \
1021 case LispComplex_t: \
1022 n##OP##_cx_cx(num, obj); \
1023 break; \
1024 default: \
1025 fatal_object_error(obj, NOT_A_NUMBER); \
1026 break; \
1027 } \
1028 } \
1029 else { \
1030 switch (OBJECT_TYPE(obj)) { \
1031 case LispFixnum_t: \
1032 switch (NRTYPE(num)) { \
1033 case N_FIXNUM: \
1034 r##OP##_fi_fi(NREAL(num), OFI(obj)); \
1035 break; \
1036 case N_BIGNUM: \
1037 r##OP##_bi_fi(NREAL(num), OFI(obj)); \
1038 break; \
1039 case N_FLONUM: \
1040 r##OP##_flonum(NREAL(num), NRFF(num), \
1041 (double)OFI(obj)); \
1042 break; \
1043 case N_FIXRATIO: \
1044 r##OP##_fr_fi(NREAL(num), OFI(obj)); \
1045 break; \
1046 case N_BIGRATIO: \
1047 r##OP##_br_fi(NREAL(num), OFI(obj)); \
1048 break; \
1049 } \
1050 break; \
1051 case LispInteger_t: \
1052 switch (NRTYPE(num)) { \
1053 case N_FIXNUM: \
1054 r##OP##_fi_fi(NREAL(num), OII(obj)); \
1055 break; \
1056 case N_BIGNUM: \
1057 r##OP##_bi_fi(NREAL(num), OII(obj)); \
1058 break; \
1059 case N_FLONUM: \
1060 r##OP##_flonum(NREAL(num), NRFF(num), \
1061 (double)OII(obj)); \
1062 break; \
1063 case N_FIXRATIO: \
1064 r##OP##_fr_fi(NREAL(num), OII(obj)); \
1065 break; \
1066 case N_BIGRATIO: \
1067 r##OP##_br_fi(NREAL(num), OII(obj)); \
1068 break; \
1069 } \
1070 break; \
1071 case LispBignum_t: \
1072 switch (NRTYPE(num)) { \
1073 case N_FIXNUM: \
1074 r##OP##_fi_bi(NREAL(num), OBI(obj)); \
1075 break; \
1076 case N_BIGNUM: \
1077 r##OP##_bi_bi(NREAL(num), OBI(obj)); \
1078 break; \
1079 case N_FLONUM: \
1080 r##OP##_flonum(NREAL(num), NRFF(num), \
1081 bi_getd(OBI(obj))); \
1082 break; \
1083 case N_FIXRATIO: \
1084 r##OP##_fr_bi(NREAL(num), OBI(obj)); \
1085 break; \
1086 case N_BIGRATIO: \
1087 r##OP##_br_bi(NREAL(num), OBI(obj)); \
1088 break; \
1089 } \
1090 break; \
1091 case LispDFloat_t: \
1092 switch (NRTYPE(num)) { \
1093 case N_FIXNUM: \
1094 r##OP##_flonum(NREAL(num), (double)NRFI(num), \
1095 ODF(obj)); \
1096 break; \
1097 case N_BIGNUM: \
1098 r##OP##_flonum(NREAL(num), bi_getd(NRBI(num)), \
1099 ODF(obj)); \
1100 break; \
1101 case N_FLONUM: \
1102 r##OP##_flonum(NREAL(num), NRFF(num), ODF(obj));\
1103 break; \
1104 case N_FIXRATIO: \
1105 r##OP##_flonum(NREAL(num), \
1106 (double)NRFRN(num) / \
1107 (double)NRFRD(num), \
1108 ODF(obj)); \
1109 break; \
1110 case N_BIGRATIO: \
1111 r##OP##_flonum(NREAL(num), br_getd(NRBR(num)), \
1112 ODF(obj)); \
1113 break; \
1114 } \
1115 break; \
1116 case LispRatio_t: \
1117 switch (NRTYPE(num)) { \
1118 case N_FIXNUM: \
1119 r##OP##_fi_fr(NREAL(num), OFRN(obj), OFRD(obj));\
1120 break; \
1121 case N_BIGNUM: \
1122 r##OP##_bi_fr(NREAL(num), OFRN(obj), OFRD(obj));\
1123 break; \
1124 case N_FLONUM: \
1125 r##OP##_flonum(NREAL(num), NRFF(num), \
1126 (double)OFRN(obj) / \
1127 (double)OFRD(obj)); \
1128 break; \
1129 case N_FIXRATIO: \
1130 r##OP##_fr_fr(NREAL(num), OFRN(obj), OFRD(obj));\
1131 break; \
1132 case N_BIGRATIO: \
1133 r##OP##_br_fr(NREAL(num), OFRN(obj), OFRD(obj));\
1134 break; \
1135 } \
1136 break; \
1137 case LispBigratio_t: \
1138 switch (NRTYPE(num)) { \
1139 case N_FIXNUM: \
1140 r##OP##_fi_br(NREAL(num), OBR(obj)); \
1141 break; \
1142 case N_BIGNUM: \
1143 r##OP##_bi_br(NREAL(num), OBR(obj)); \
1144 break; \
1145 case N_FLONUM: \
1146 r##OP##_flonum(NREAL(num), NRFF(num), \
1147 br_getd(OBR(obj))); \
1148 break; \
1149 case N_FIXRATIO: \
1150 r##OP##_fr_br(NREAL(num), OBR(obj)); \
1151 break; \
1152 case N_BIGRATIO: \
1153 r##OP##_br_br(NREAL(num), OBR(obj)); \
1154 break; \
1155 } \
1156 break; \
1157 case LispComplex_t: \
1158 n##OP##_re_cx(num, obj); \
1159 break; \
1160 default: \
1161 fatal_object_error(obj, NOT_A_NUMBER); \
1162 break; \
1163 } \
1164 } \
1165 }
1166
1167 static void
1168 DEFOP_NUMBER_OBJECT(add)
1169
1170 static void
1171 DEFOP_NUMBER_OBJECT(sub)
1172
1173 static void
1174 DEFOP_NUMBER_OBJECT(div)
1175
1176 static void
1177 DEFOP_NUMBER_OBJECT(mul)
1178
1179
1180 /************************************************************************
1181 * ABS
1182 ************************************************************************/
1183 static void
1184 abs_real(n_real *real)
1185 {
1186 switch (RTYPE(real)) {
1187 case N_FIXNUM: rabs_fi(real); break;
1188 case N_BIGNUM: rabs_bi(real); break;
1189 case N_FLONUM: rabs_ff(real); break;
1190 case N_FIXRATIO: rabs_fr(real); break;
1191 case N_BIGRATIO: rabs_br(real); break;
1192 }
1193 }
1194
1195 static void
abs_number(n_number * num)1196 abs_number(n_number *num)
1197 {
1198 if (num->complex)
1199 nabs_cx(num);
1200 else {
1201 switch (NRTYPE(num)) {
1202 case N_FIXNUM: rabs_fi(NREAL(num)); break;
1203 case N_BIGNUM: rabs_bi(NREAL(num)); break;
1204 case N_FLONUM: rabs_ff(NREAL(num)); break;
1205 case N_FIXRATIO: rabs_fr(NREAL(num)); break;
1206 case N_BIGRATIO: rabs_br(NREAL(num)); break;
1207 }
1208 }
1209 }
1210
1211 static void
nabs_cx(n_number * num)1212 nabs_cx(n_number *num)
1213 {
1214 n_real temp;
1215
1216 abs_real(NREAL(num));
1217 abs_real(NIMAG(num));
1218
1219 if (cmp_real_real(NREAL(num), NIMAG(num)) < 0) {
1220 memcpy(&temp, NIMAG(num), sizeof(n_real));
1221 memcpy(NIMAG(num), NREAL(num), sizeof(n_real));
1222 memcpy(NREAL(num), &temp, sizeof(n_real));
1223 }
1224
1225 if (cmp_real_real(NIMAG(num), &zero) == 0) {
1226 num->complex = 0;
1227 if (NITYPE(num) == N_FLONUM) {
1228 /* change number type */
1229 temp.type = N_FLONUM;
1230 temp.data.flonum = 1.0;
1231 mul_real_real(NREAL(num), &temp);
1232 }
1233 else
1234 clear_real(NIMAG(num));
1235 }
1236 else {
1237 div_real_real(NIMAG(num), NREAL(num));
1238 set_real_real(&temp, NIMAG(num));
1239 mul_real_real(NIMAG(num), &temp);
1240 clear_real(&temp);
1241
1242 add_real_real(NIMAG(num), &one);
1243 sqrt_real(NIMAG(num));
1244
1245 mul_real_real(NIMAG(num), NREAL(num));
1246 clear_real(NREAL(num));
1247 memcpy(NREAL(num), NIMAG(num), sizeof(n_real));
1248 num->complex = 0;
1249 }
1250 }
1251
1252 static INLINE void
rabs_fi(n_real * real)1253 rabs_fi(n_real *real)
1254 {
1255 if (RFI(real) < 0)
1256 rneg_fi(real);
1257 }
1258
1259 static INLINE void
rabs_bi(n_real * real)1260 rabs_bi(n_real *real)
1261 {
1262 if (mpi_cmpi(RBI(real), 0) < 0)
1263 mpi_neg(RBI(real), RBI(real));
1264 }
1265
1266 static INLINE void
rabs_ff(n_real * real)1267 rabs_ff(n_real *real)
1268 {
1269 if (RFF(real) < 0.0)
1270 RFF(real) = -RFF(real);
1271 }
1272
1273 static INLINE void
rabs_fr(n_real * real)1274 rabs_fr(n_real *real)
1275 {
1276 if (RFRN(real) < 0)
1277 rneg_fr(real);
1278 }
1279
1280 static INLINE void
rabs_br(n_real * real)1281 rabs_br(n_real *real)
1282 {
1283 if (mpi_cmpi(RBRN(real), 0) < 0)
1284 mpi_neg(RBRN(real), RBRN(real));
1285 }
1286
1287
1288 /************************************************************************
1289 * NEG
1290 ************************************************************************/
1291 static void
neg_real(n_real * real)1292 neg_real(n_real *real)
1293 {
1294 switch (RTYPE(real)) {
1295 case N_FIXNUM: rneg_fi(real); break;
1296 case N_BIGNUM: rneg_bi(real); break;
1297 case N_FLONUM: rneg_ff(real); break;
1298 case N_FIXRATIO: rneg_fr(real); break;
1299 case N_BIGRATIO: rneg_br(real); break;
1300 }
1301 }
1302
1303 static void
neg_number(n_number * num)1304 neg_number(n_number *num)
1305 {
1306 if (num->complex) {
1307 neg_real(NREAL(num));
1308 neg_real(NIMAG(num));
1309 }
1310 else {
1311 switch (NRTYPE(num)) {
1312 case N_FIXNUM: rneg_fi(NREAL(num)); break;
1313 case N_BIGNUM: rneg_bi(NREAL(num)); break;
1314 case N_FLONUM: rneg_ff(NREAL(num)); break;
1315 case N_FIXRATIO: rneg_fr(NREAL(num)); break;
1316 case N_BIGRATIO: rneg_br(NREAL(num)); break;
1317 }
1318 }
1319 }
1320
1321 static void
rneg_fi(n_real * real)1322 rneg_fi(n_real *real)
1323 {
1324 if (RFI(real) == MINSLONG) {
1325 mpi *bigi = XALLOC(mpi);
1326
1327 mpi_init(bigi);
1328 mpi_seti(bigi, RFI(real));
1329 mpi_neg(bigi, bigi);
1330 RTYPE(real) = N_BIGNUM;
1331 RBI(real) = bigi;
1332 }
1333 else
1334 RFI(real) = -RFI(real);
1335 }
1336
1337 static INLINE void
rneg_bi(n_real * real)1338 rneg_bi(n_real *real)
1339 {
1340 mpi_neg(RBI(real), RBI(real));
1341 }
1342
1343 static INLINE void
rneg_ff(n_real * real)1344 rneg_ff(n_real *real)
1345 {
1346 RFF(real) = -RFF(real);
1347 }
1348
1349 static void
rneg_fr(n_real * real)1350 rneg_fr(n_real *real)
1351 {
1352 if (RFRN(real) == MINSLONG) {
1353 mpr *bigr = XALLOC(mpr);
1354
1355 mpr_init(bigr);
1356 mpr_seti(bigr, RFRN(real), RFRD(real));
1357 mpi_neg(mpr_num(bigr), mpr_num(bigr));
1358 RTYPE(real) = N_BIGRATIO;
1359 RBR(real) = bigr;
1360 }
1361 else
1362 RFRN(real) = -RFRN(real);
1363 }
1364
1365 static INLINE void
rneg_br(n_real * real)1366 rneg_br(n_real *real)
1367 {
1368 mpi_neg(RBRN(real), RBRN(real));
1369 }
1370
1371
1372 /************************************************************************
1373 * SQRT
1374 ************************************************************************/
1375 static void
sqrt_real(n_real * real)1376 sqrt_real(n_real *real)
1377 {
1378 switch (RTYPE(real)) {
1379 case N_FIXNUM:
1380 case N_BIGNUM:
1381 rsqrt_xi(real);
1382 break;
1383 case N_FLONUM:
1384 rsqrt_ff(real);
1385 break;
1386 case N_FIXRATIO:
1387 case N_BIGRATIO:
1388 rsqrt_xr(real);
1389 break;
1390 }
1391 }
1392
1393 static void
sqrt_number(n_number * num)1394 sqrt_number(n_number *num)
1395 {
1396 if (num->complex)
1397 nsqrt_cx(num);
1398 else {
1399 switch (NRTYPE(num)) {
1400 case N_FIXNUM:
1401 case N_BIGNUM:
1402 nsqrt_xi(num);
1403 break;
1404 case N_FLONUM:
1405 nsqrt_ff(num);
1406 break;
1407 case N_FIXRATIO:
1408 case N_BIGRATIO:
1409 nsqrt_xr(num);
1410 break;
1411 }
1412 }
1413 }
1414
1415 static void
rsqrt_xi(n_real * real)1416 rsqrt_xi(n_real *real)
1417 {
1418 int exact;
1419 mpi bignum;
1420
1421 if (cmp_real_real(real, &zero) < 0)
1422 fatal_error(FLOATING_POINT_EXCEPTION);
1423
1424 mpi_init(&bignum);
1425 if (RTYPE(real) == N_BIGNUM)
1426 exact = mpi_sqrt(&bignum, RBI(real));
1427 else {
1428 mpi tmp;
1429
1430 mpi_init(&tmp);
1431 mpi_seti(&tmp, RFI(real));
1432 exact = mpi_sqrt(&bignum, &tmp);
1433 mpi_clear(&tmp);
1434 }
1435 if (exact) {
1436 if (RTYPE(real) == N_BIGNUM) {
1437 mpi_set(RBI(real), &bignum);
1438 rbi_canonicalize(real);
1439 }
1440 else
1441 RFI(real) = mpi_geti(&bignum);
1442 }
1443 else {
1444 double value;
1445
1446 if (RTYPE(real) == N_BIGNUM) {
1447 value = bi_getd(RBI(real));
1448 RCLEAR_BI(real);
1449 }
1450 else
1451 value = (double)RFI(real);
1452
1453 value = sqrt(value);
1454 RTYPE(real) = N_FLONUM;
1455 RFF(real) = value;
1456 }
1457 mpi_clear(&bignum);
1458 }
1459
1460 static void
rsqrt_xr(n_real * real)1461 rsqrt_xr(n_real *real)
1462 {
1463 n_real num, den;
1464
1465 if (cmp_real_real(real, &zero) < 0)
1466 fatal_error(FLOATING_POINT_EXCEPTION);
1467
1468 if (RTYPE(real) == N_FIXRATIO) {
1469 num.type = den.type = N_FIXNUM;
1470 num.data.fixnum = RFRN(real);
1471 den.data.fixnum = RFRD(real);
1472 }
1473 else {
1474 mpi *bignum;
1475
1476 if (mpi_fiti(RBRN(real))) {
1477 num.type = N_FIXNUM;
1478 num.data.fixnum = mpi_geti(RBRN(real));
1479 }
1480 else {
1481 bignum = XALLOC(mpi);
1482 mpi_init(bignum);
1483 mpi_set(bignum, RBRN(real));
1484 num.type = N_BIGNUM;
1485 num.data.bignum = bignum;
1486 }
1487
1488 if (mpi_fiti(RBRD(real))) {
1489 den.type = N_FIXNUM;
1490 den.data.fixnum = mpi_geti(RBRD(real));
1491 }
1492 else {
1493 bignum = XALLOC(mpi);
1494 mpi_init(bignum);
1495 mpi_set(bignum, RBRD(real));
1496 den.type = N_BIGNUM;
1497 den.data.bignum = bignum;
1498 }
1499 }
1500
1501 rsqrt_xi(&num);
1502 rsqrt_xi(&den);
1503
1504 clear_real(real);
1505 memcpy(real, &num, sizeof(n_real));
1506 div_real_real(real, &den);
1507 clear_real(&den);
1508 }
1509
1510 static void
rsqrt_ff(n_real * real)1511 rsqrt_ff(n_real *real)
1512 {
1513 if (RFF(real) < 0.0)
1514 fatal_error(FLOATING_POINT_EXCEPTION);
1515 RFF(real) = sqrt(RFF(real));
1516 }
1517
1518
1519 static void
nsqrt_cx(n_number * num)1520 nsqrt_cx(n_number *num)
1521 {
1522 n_number mag;
1523 n_real *real, *imag;
1524
1525 real = &(mag.real);
1526 imag = &(mag.imag);
1527 set_real_real(real, NREAL(num));
1528 set_real_real(imag, NIMAG(num));
1529 mag.complex = 1;
1530
1531 nabs_cx(&mag); /* this will free the imag part data */
1532 if (cmp_real_real(real, &zero) == 0) {
1533 clear_number(num);
1534 memcpy(NREAL(num), real, sizeof(n_real));
1535 clear_real(real);
1536 num->complex = 0;
1537 return;
1538 }
1539 else if (cmp_real_real(NREAL(num), &zero) > 0) {
1540 /* R = sqrt((mag + Ra) / 2) */
1541 add_real_real(NREAL(num), real);
1542 clear_real(real);
1543 div_real_real(NREAL(num), &two);
1544 sqrt_real(NREAL(num));
1545
1546 /* I = Ia / R / 2 */
1547 div_real_real(NIMAG(num), NREAL(num));
1548 div_real_real(NIMAG(num), &two);
1549 }
1550 else {
1551 /* remember old imag part */
1552 memcpy(imag, NIMAG(num), sizeof(n_real));
1553
1554 /* I = sqrt((mag - Ra) / 2) */
1555 memcpy(NIMAG(num), real, sizeof(n_real));
1556 sub_real_real(NIMAG(num), NREAL(num));
1557 div_real_real(NIMAG(num), &two);
1558 sqrt_real(NIMAG(num));
1559 if (cmp_real_real(imag, &zero) < 0)
1560 neg_real(NIMAG(num));
1561
1562 /* R = Ia / I / 2 */
1563 clear_real(NREAL(num));
1564 /* start with old imag part */
1565 memcpy(NREAL(num), imag, sizeof(n_real));
1566 div_real_real(NREAL(num), NIMAG(num));
1567 div_real_real(NREAL(num), &two);
1568 }
1569
1570 ncx_canonicalize(num);
1571 }
1572
1573 static void
nsqrt_xi(n_number * num)1574 nsqrt_xi(n_number *num)
1575 {
1576 if (cmp_real_real(NREAL(num), &zero) < 0) {
1577 memcpy(NIMAG(num), NREAL(num), sizeof(n_real));
1578 neg_real(NIMAG(num));
1579 rsqrt_xi(NIMAG(num));
1580 NRTYPE(num) = N_FIXNUM;
1581 NRFI(num) = 0;
1582 num->complex = 1;
1583 }
1584 else
1585 rsqrt_xi(NREAL(num));
1586 }
1587
1588 static void
nsqrt_ff(n_number * num)1589 nsqrt_ff(n_number *num)
1590 {
1591 double value;
1592
1593 if (NRFF(num) < 0.0) {
1594 value = sqrt(-NRFF(num));
1595
1596 NITYPE(num) = N_FLONUM;
1597 NIFF(num) = value;
1598 NRTYPE(num) = N_FIXNUM;
1599 NRFI(num) = 0;
1600 num->complex = 1;
1601 }
1602 else {
1603 value = sqrt(NRFF(num));
1604 NRFF(num) = value;
1605 }
1606 }
1607
1608 static void
nsqrt_xr(n_number * num)1609 nsqrt_xr(n_number *num)
1610 {
1611 if (cmp_real_real(NREAL(num), &zero) < 0) {
1612 memcpy(NIMAG(num), NREAL(num), sizeof(n_real));
1613 neg_real(NIMAG(num));
1614 rsqrt_xr(NIMAG(num));
1615 NRTYPE(num) = N_FIXNUM;
1616 NRFI(num) = 0;
1617 num->complex = 1;
1618 }
1619 else
1620 rsqrt_xr(NREAL(num));
1621 }
1622
1623
1624 /************************************************************************
1625 * MOD
1626 ************************************************************************/
1627 static void
mod_real_real(n_real * real,n_real * val)1628 mod_real_real(n_real *real, n_real *val)
1629 {
1630 /* Assume both operands are integers */
1631 switch (RTYPE(real)) {
1632 case N_FIXNUM:
1633 switch (RTYPE(val)) {
1634 case N_FIXNUM:
1635 rmod_fi_fi(real, RFI(val));
1636 break;
1637 case N_BIGNUM:
1638 rmod_fi_bi(real, RBI(val));
1639 break;
1640 }
1641 break;
1642 case N_BIGNUM:
1643 switch (RTYPE(val)) {
1644 case N_FIXNUM:
1645 rmod_bi_fi(real, RFI(val));
1646 break;
1647 case N_BIGNUM:
1648 rmod_bi_bi(real, RBI(val));
1649 break;
1650 }
1651 break;
1652 }
1653 }
1654
1655 static void
mod_real_object(n_real * real,LispObj * obj)1656 mod_real_object(n_real *real, LispObj *obj)
1657 {
1658 switch (RTYPE(real)) {
1659 case N_FIXNUM:
1660 switch (OBJECT_TYPE(obj)) {
1661 case LispFixnum_t:
1662 rmod_fi_fi(real, OFI(obj));
1663 return;
1664 case LispInteger_t:
1665 rmod_fi_fi(real, OII(obj));
1666 return;
1667 case LispBignum_t:
1668 rmod_fi_bi(real, OBI(obj));
1669 return;
1670 default:
1671 break;
1672 }
1673 break;
1674 case N_BIGNUM:
1675 switch (OBJECT_TYPE(obj)) {
1676 case LispFixnum_t:
1677 rmod_bi_fi(real, OFI(obj));
1678 return;
1679 case LispInteger_t:
1680 rmod_bi_fi(real, OII(obj));
1681 return;
1682 case LispBignum_t:
1683 rmod_bi_bi(real, OBI(obj));
1684 return;
1685 default:
1686 break;
1687 }
1688 break;
1689 /* Assume the n_real object is an integer */
1690 }
1691 fatal_object_error(obj, NOT_AN_INTEGER);
1692 }
1693
1694 static void
rmod_fi_fi(n_real * real,long fi)1695 rmod_fi_fi(n_real *real, long fi)
1696 {
1697 if (fi == 0)
1698 fatal_error(DIVIDE_BY_ZERO);
1699
1700 if ((RFI(real) < 0) ^ (fi < 0))
1701 RFI(real) = (RFI(real) % fi) + fi;
1702 else if (RFI(real) == MINSLONG || fi == MINSLONG) {
1703 mpi bignum;
1704
1705 mpi_init(&bignum);
1706 mpi_seti(&bignum, RFI(real));
1707 RFI(real) = mpi_modi(&bignum, fi);
1708 mpi_clear(&bignum);
1709 }
1710 else
1711 RFI(real) = RFI(real) % fi;
1712 }
1713
1714 static void
rmod_fi_bi(n_real * real,mpi * bignum)1715 rmod_fi_bi(n_real *real, mpi *bignum)
1716 {
1717 mpi *bigi;
1718
1719 if (mpi_cmpi(bignum, 0) == 0)
1720 fatal_error(DIVIDE_BY_ZERO);
1721
1722 bigi = XALLOC(mpi);
1723 mpi_init(bigi);
1724 mpi_seti(bigi, RFI(real));
1725 mpi_mod(bigi, bigi, bignum);
1726 RTYPE(real) = N_BIGNUM;
1727 RBI(real) = bigi;
1728 rbi_canonicalize(real);
1729 }
1730
1731 static void
rmod_bi_fi(n_real * real,long fi)1732 rmod_bi_fi(n_real *real, long fi)
1733 {
1734 mpi iop;
1735
1736 if (fi == 0)
1737 fatal_error(DIVIDE_BY_ZERO);
1738
1739 mpi_init(&iop);
1740 mpi_seti(&iop, fi);
1741 mpi_mod(RBI(real), RBI(real), &iop);
1742 mpi_clear(&iop);
1743 rbi_canonicalize(real);
1744 }
1745
1746 static void
rmod_bi_bi(n_real * real,mpi * bignum)1747 rmod_bi_bi(n_real *real, mpi *bignum)
1748 {
1749 if (mpi_cmpi(bignum, 0) == 0)
1750 fatal_error(DIVIDE_BY_ZERO);
1751
1752 mpi_mod(RBI(real), RBI(real), bignum);
1753 rbi_canonicalize(real);
1754 }
1755
1756 /************************************************************************
1757 * REM
1758 ************************************************************************/
1759 static void
rem_real_object(n_real * real,LispObj * obj)1760 rem_real_object(n_real *real, LispObj *obj)
1761 {
1762 switch (RTYPE(real)) {
1763 case N_FIXNUM:
1764 switch (OBJECT_TYPE(obj)) {
1765 case LispFixnum_t:
1766 rrem_fi_fi(real, OFI(obj));
1767 return;
1768 case LispInteger_t:
1769 rrem_fi_fi(real, OII(obj));
1770 return;
1771 case LispBignum_t:
1772 rrem_fi_bi(real, OBI(obj));
1773 return;
1774 default:
1775 break;
1776 }
1777 break;
1778 case N_BIGNUM:
1779 switch (OBJECT_TYPE(obj)) {
1780 case LispFixnum_t:
1781 rrem_bi_fi(real, OFI(obj));
1782 return;
1783 case LispInteger_t:
1784 rrem_bi_fi(real, OII(obj));
1785 return;
1786 case LispBignum_t:
1787 rrem_bi_bi(real, OBI(obj));
1788 return;
1789 default:
1790 break;
1791 }
1792 break;
1793 /* Assume the n_real object is an integer */
1794 }
1795 fatal_object_error(obj, NOT_AN_INTEGER);
1796 }
1797
1798 static void
rrem_fi_fi(n_real * real,long fi)1799 rrem_fi_fi(n_real *real, long fi)
1800 {
1801 if (fi == 0)
1802 fatal_error(DIVIDE_BY_ZERO);
1803
1804 if (RFI(real) == MINSLONG || fi == MINSLONG) {
1805 mpi bignum;
1806
1807 mpi_init(&bignum);
1808 mpi_seti(&bignum, RFI(real));
1809 RFI(real) = mpi_remi(&bignum, fi);
1810 mpi_clear(&bignum);
1811 }
1812 else
1813 RFI(real) = RFI(real) % fi;
1814 }
1815
1816 static void
rrem_fi_bi(n_real * real,mpi * bignum)1817 rrem_fi_bi(n_real *real, mpi *bignum)
1818 {
1819 mpi *bigi;
1820
1821 if (mpi_cmpi(bignum, 0) == 0)
1822 fatal_error(DIVIDE_BY_ZERO);
1823
1824 bigi = XALLOC(mpi);
1825 mpi_init(bigi);
1826 mpi_seti(bigi, RFI(real));
1827 mpi_rem(bigi, bigi, bignum);
1828 RTYPE(real) = N_BIGNUM;
1829 RBI(real) = bigi;
1830 rbi_canonicalize(real);
1831 }
1832
1833 static void
rrem_bi_fi(n_real * real,long fi)1834 rrem_bi_fi(n_real *real, long fi)
1835 {
1836 mpi iop;
1837
1838 if (fi == 0)
1839 fatal_error(DIVIDE_BY_ZERO);
1840
1841 mpi_init(&iop);
1842 mpi_seti(&iop, fi);
1843 mpi_rem(RBI(real), RBI(real), &iop);
1844 mpi_clear(&iop);
1845 rbi_canonicalize(real);
1846 }
1847
1848 static void
rrem_bi_bi(n_real * real,mpi * bignum)1849 rrem_bi_bi(n_real *real, mpi *bignum)
1850 {
1851 if (mpi_cmpi(bignum, 0) == 0)
1852 fatal_error(DIVIDE_BY_ZERO);
1853
1854 mpi_rem(RBI(real), RBI(real), bignum);
1855 rbi_canonicalize(real);
1856 }
1857
1858
1859 /************************************************************************
1860 * GCD
1861 ************************************************************************/
1862 static void
gcd_real_object(n_real * real,LispObj * obj)1863 gcd_real_object(n_real *real, LispObj *obj)
1864 {
1865 if (!INTEGERP(obj))
1866 fatal_object_error(obj, NOT_AN_INTEGER);
1867
1868 /* check for zero operand */
1869 if (cmp_real_real(real, &zero) == 0)
1870 set_real_object(real, obj);
1871 else if (cmp_real_object(&zero, obj) != 0) {
1872 n_real rest, temp;
1873
1874 set_real_object(&rest, obj);
1875 for (;;) {
1876 mod_real_real(&rest, real);
1877 if (cmp_real_real(&rest, &zero) == 0)
1878 break;
1879 memcpy(&temp, real, sizeof(n_real));
1880 memcpy(real, &rest, sizeof(n_real));
1881 memcpy(&rest, &temp, sizeof(n_real));
1882 }
1883 clear_real(&rest);
1884 }
1885 }
1886
1887 /************************************************************************
1888 * AND
1889 ************************************************************************/
1890 static void
and_real_object(n_real * real,LispObj * obj)1891 and_real_object(n_real *real, LispObj *obj)
1892 {
1893 mpi *bigi, iop;
1894
1895 switch (OBJECT_TYPE(obj)) {
1896 case LispFixnum_t:
1897 switch (RTYPE(real)) {
1898 case N_FIXNUM:
1899 RFI(real) &= OFI(obj);
1900 break;
1901 case N_BIGNUM:
1902 mpi_init(&iop);
1903 mpi_seti(&iop, OFI(obj));
1904 mpi_and(RBI(real), RBI(real), &iop);
1905 mpi_clear(&iop);
1906 rbi_canonicalize(real);
1907 break;
1908 }
1909 break;
1910 case LispInteger_t:
1911 switch (RTYPE(real)) {
1912 case N_FIXNUM:
1913 RFI(real) &= OII(obj);
1914 break;
1915 case N_BIGNUM:
1916 mpi_init(&iop);
1917 mpi_seti(&iop, OII(obj));
1918 mpi_and(RBI(real), RBI(real), &iop);
1919 mpi_clear(&iop);
1920 rbi_canonicalize(real);
1921 break;
1922 }
1923 break;
1924 case LispBignum_t:
1925 switch (RTYPE(real)) {
1926 case N_FIXNUM:
1927 bigi = XALLOC(mpi);
1928 mpi_init(bigi);
1929 mpi_seti(bigi, RFI(real));
1930 mpi_and(bigi, bigi, OBI(obj));
1931 RTYPE(real) = N_BIGNUM;
1932 RBI(real) = bigi;
1933 rbi_canonicalize(real);
1934 break;
1935 case N_BIGNUM:
1936 mpi_and(RBI(real), RBI(real), OBI(obj));
1937 rbi_canonicalize(real);
1938 break;
1939 }
1940 break;
1941 default:
1942 fatal_object_error(obj, NOT_AN_INTEGER);
1943 break;
1944 }
1945 }
1946
1947
1948 /************************************************************************
1949 * EQV
1950 ************************************************************************/
1951 static void
eqv_real_object(n_real * real,LispObj * obj)1952 eqv_real_object(n_real *real, LispObj *obj)
1953 {
1954 mpi *bigi, iop;
1955
1956 switch (OBJECT_TYPE(obj)) {
1957 case LispFixnum_t:
1958 switch (RTYPE(real)) {
1959 case N_FIXNUM:
1960 RFI(real) ^= ~OFI(obj);
1961 break;
1962 case N_BIGNUM:
1963 mpi_init(&iop);
1964 mpi_seti(&iop, OFI(obj));
1965 mpi_com(&iop, &iop);
1966 mpi_xor(RBI(real), RBI(real), &iop);
1967 mpi_clear(&iop);
1968 rbi_canonicalize(real);
1969 break;
1970 }
1971 break;
1972 case LispInteger_t:
1973 switch (RTYPE(real)) {
1974 case N_FIXNUM:
1975 RFI(real) ^= ~OII(obj);
1976 break;
1977 case N_BIGNUM:
1978 mpi_init(&iop);
1979 mpi_seti(&iop, OII(obj));
1980 mpi_com(&iop, &iop);
1981 mpi_xor(RBI(real), RBI(real), &iop);
1982 mpi_clear(&iop);
1983 rbi_canonicalize(real);
1984 break;
1985 }
1986 break;
1987 case LispBignum_t:
1988 switch (RTYPE(real)) {
1989 case N_FIXNUM:
1990 bigi = XALLOC(mpi);
1991 mpi_init(bigi);
1992 mpi_seti(bigi, RFI(real));
1993 mpi_com(bigi, bigi);
1994 mpi_xor(bigi, bigi, OBI(obj));
1995 RTYPE(real) = N_BIGNUM;
1996 RBI(real) = bigi;
1997 rbi_canonicalize(real);
1998 break;
1999 case N_BIGNUM:
2000 mpi_com(RBI(real), RBI(real));
2001 mpi_xor(RBI(real), RBI(real), OBI(obj));
2002 rbi_canonicalize(real);
2003 break;
2004 }
2005 break;
2006 default:
2007 fatal_object_error(obj, NOT_AN_INTEGER);
2008 break;
2009 }
2010 }
2011
2012
2013 /************************************************************************
2014 * IOR
2015 ************************************************************************/
2016 static void
ior_real_object(n_real * real,LispObj * obj)2017 ior_real_object(n_real *real, LispObj *obj)
2018 {
2019 mpi *bigi, iop;
2020
2021 switch (OBJECT_TYPE(obj)) {
2022 case LispFixnum_t:
2023 switch (RTYPE(real)) {
2024 case N_FIXNUM:
2025 RFI(real) |= OFI(obj);
2026 break;
2027 case N_BIGNUM:
2028 mpi_init(&iop);
2029 mpi_seti(&iop, OFI(obj));
2030 mpi_ior(RBI(real), RBI(real), &iop);
2031 mpi_clear(&iop);
2032 rbi_canonicalize(real);
2033 break;
2034 }
2035 break;
2036 case LispInteger_t:
2037 switch (RTYPE(real)) {
2038 case N_FIXNUM:
2039 RFI(real) |= OII(obj);
2040 break;
2041 case N_BIGNUM:
2042 mpi_init(&iop);
2043 mpi_seti(&iop, OII(obj));
2044 mpi_ior(RBI(real), RBI(real), &iop);
2045 mpi_clear(&iop);
2046 rbi_canonicalize(real);
2047 break;
2048 }
2049 break;
2050 case LispBignum_t:
2051 switch (RTYPE(real)) {
2052 case N_FIXNUM:
2053 bigi = XALLOC(mpi);
2054 mpi_init(bigi);
2055 mpi_seti(bigi, RFI(real));
2056 mpi_ior(bigi, bigi, OBI(obj));
2057 RTYPE(real) = N_BIGNUM;
2058 RBI(real) = bigi;
2059 rbi_canonicalize(real);
2060 break;
2061 case N_BIGNUM:
2062 mpi_ior(RBI(real), RBI(real), OBI(obj));
2063 rbi_canonicalize(real);
2064 break;
2065 }
2066 break;
2067 default:
2068 fatal_object_error(obj, NOT_AN_INTEGER);
2069 break;
2070 }
2071 }
2072
2073
2074 /************************************************************************
2075 * NOT
2076 ************************************************************************/
2077 static void
not_real(n_real * real)2078 not_real(n_real *real)
2079 {
2080 if (RTYPE(real) == N_FIXNUM)
2081 RFI(real) = ~RFI(real);
2082 else {
2083 mpi_com(RBI(real), RBI(real));
2084 rbi_canonicalize(real);
2085 }
2086 }
2087
2088 /************************************************************************
2089 * XOR
2090 ************************************************************************/
2091 static void
xor_real_object(n_real * real,LispObj * obj)2092 xor_real_object(n_real *real, LispObj *obj)
2093 {
2094 mpi *bigi, iop;
2095
2096 switch (OBJECT_TYPE(obj)) {
2097 case LispFixnum_t:
2098 switch (RTYPE(real)) {
2099 case N_FIXNUM:
2100 RFI(real) ^= OFI(obj);
2101 break;
2102 case N_BIGNUM:
2103 mpi_init(&iop);
2104 mpi_seti(&iop, OFI(obj));
2105 mpi_xor(RBI(real), RBI(real), &iop);
2106 mpi_clear(&iop);
2107 rbi_canonicalize(real);
2108 break;
2109 }
2110 break;
2111 case LispInteger_t:
2112 switch (RTYPE(real)) {
2113 case N_FIXNUM:
2114 RFI(real) ^= OII(obj);
2115 break;
2116 case N_BIGNUM:
2117 mpi_init(&iop);
2118 mpi_seti(&iop, OII(obj));
2119 mpi_xor(RBI(real), RBI(real), &iop);
2120 mpi_clear(&iop);
2121 rbi_canonicalize(real);
2122 break;
2123 }
2124 break;
2125 case LispBignum_t:
2126 switch (RTYPE(real)) {
2127 case N_FIXNUM:
2128 bigi = XALLOC(mpi);
2129 mpi_init(bigi);
2130 mpi_seti(bigi, RFI(real));
2131 mpi_xor(bigi, bigi, OBI(obj));
2132 RTYPE(real) = N_BIGNUM;
2133 RBI(real) = bigi;
2134 rbi_canonicalize(real);
2135 break;
2136 case N_BIGNUM:
2137 mpi_xor(RBI(real), RBI(real), OBI(obj));
2138 rbi_canonicalize(real);
2139 break;
2140 }
2141 break;
2142 default:
2143 fatal_object_error(obj, NOT_AN_INTEGER);
2144 break;
2145 }
2146 }
2147
2148
2149 /************************************************************************
2150 * DIVIDE
2151 ************************************************************************/
2152 static void
divide_number_object(n_number * num,LispObj * obj,int fun,int flo)2153 divide_number_object(n_number *num, LispObj *obj, int fun, int flo)
2154 {
2155 switch (OBJECT_TYPE(obj)) {
2156 case LispFixnum_t:
2157 switch (NRTYPE(num)) {
2158 case N_FIXNUM:
2159 ndivide_fi_fi(num, OFI(obj), fun, flo);
2160 break;
2161 case N_BIGNUM:
2162 ndivide_xi_xi(num, obj, fun, flo);
2163 break;
2164 case N_FLONUM:
2165 ndivide_flonum(num, NRFF(num), (double)OFI(obj), fun, flo);
2166 break;
2167 case N_FIXRATIO:
2168 case N_BIGRATIO:
2169 ndivide_xr_xi(num, obj, fun, flo);
2170 break;
2171 }
2172 break;
2173 case LispInteger_t:
2174 switch (NRTYPE(num)) {
2175 case N_FIXNUM:
2176 ndivide_fi_fi(num, OII(obj), fun, flo);
2177 break;
2178 case N_BIGNUM:
2179 ndivide_xi_xi(num, obj, fun, flo);
2180 break;
2181 case N_FLONUM:
2182 ndivide_flonum(num, NRFF(num), (double)OII(obj), fun, flo);
2183 break;
2184 case N_FIXRATIO:
2185 case N_BIGRATIO:
2186 ndivide_xr_xi(num, obj, fun, flo);
2187 break;
2188 }
2189 break;
2190 case LispBignum_t:
2191 switch (NRTYPE(num)) {
2192 case N_FIXNUM:
2193 case N_BIGNUM:
2194 ndivide_xi_xi(num, obj, fun, flo);
2195 break;
2196 case N_FLONUM:
2197 ndivide_flonum(num, NRFF(num), bi_getd(OBI(obj)),
2198 fun, flo);
2199 break;
2200 case N_FIXRATIO:
2201 case N_BIGRATIO:
2202 ndivide_xr_xi(num, obj, fun, flo);
2203 break;
2204 }
2205 break;
2206 case LispDFloat_t:
2207 switch (NRTYPE(num)) {
2208 case N_FIXNUM:
2209 ndivide_flonum(num, (double)NRFI(num), ODF(obj),
2210 fun, flo);
2211 break;
2212 case N_BIGNUM:
2213 ndivide_flonum(num, bi_getd(NRBI(num)), ODF(obj),
2214 fun, flo);
2215 break;
2216 case N_FLONUM:
2217 ndivide_flonum(num, NRFF(num), ODF(obj), fun, flo);
2218 break;
2219 case N_FIXRATIO:
2220 ndivide_flonum(num,
2221 (double)NRFRN(num) / (double)NRFRD(num),
2222 ODF(obj), fun, flo);
2223 break;
2224 case N_BIGRATIO:
2225 ndivide_flonum(num, br_getd(NRBR(num)), ODF(obj),
2226 fun, flo);
2227 break;
2228 }
2229 break;
2230 case LispRatio_t:
2231 switch (NRTYPE(num)) {
2232 case N_FIXNUM:
2233 case N_BIGNUM:
2234 ndivide_xi_xr(num, obj, fun, flo);
2235 break;
2236 case N_FLONUM:
2237 ndivide_flonum(num, NRFF(num),
2238 (double)OFRN(obj) / (double)OFRD(obj),
2239 fun, flo);
2240 break;
2241 case N_FIXRATIO:
2242 case N_BIGRATIO:
2243 ndivide_xr_xr(num, obj, fun, flo);
2244 break;
2245 }
2246 break;
2247 case LispBigratio_t:
2248 switch (NRTYPE(num)) {
2249 case N_FIXNUM:
2250 case N_BIGNUM:
2251 ndivide_xi_xr(num, obj, fun, flo);
2252 break;
2253 case N_FLONUM:
2254 ndivide_flonum(num, NRFF(num), br_getd(OBR(obj)),
2255 fun, flo);
2256 break;
2257 case N_FIXRATIO:
2258 case N_BIGRATIO:
2259 ndivide_xr_xr(num, obj, fun, flo);
2260 break;
2261 }
2262 break;
2263 default:
2264 fatal_object_error(obj, NOT_A_REAL_NUMBER);
2265 break;
2266 }
2267 }
2268
2269
2270 /************************************************************************
2271 * COMPARE
2272 ************************************************************************/
2273 static int
cmp_real_real(n_real * op1,n_real * op2)2274 cmp_real_real(n_real *op1, n_real *op2)
2275 {
2276 switch (RTYPE(op1)) {
2277 case N_FIXNUM:
2278 switch (RTYPE(op2)) {
2279 case N_FIXNUM:
2280 return (cmp_fi_fi(RFI(op1), RFI(op2)));
2281 case N_BIGNUM:
2282 return (cmp_fi_bi(RFI(op1), RBI(op2)));
2283 case N_FLONUM:
2284 return (cmp_flonum((double)RFI(op1), RFF(op2)));
2285 case N_FIXRATIO:
2286 return (cmp_fi_fr(RFI(op1), RFRN(op2), RFRD(op2)));
2287 case N_BIGRATIO:
2288 return (cmp_fi_br(RFI(op1), RBR(op2)));
2289 }
2290 break;
2291 case N_BIGNUM:
2292 switch (RTYPE(op2)) {
2293 case N_FIXNUM:
2294 return (cmp_bi_fi(RBI(op1), RFI(op2)));
2295 case N_BIGNUM:
2296 return (cmp_bi_bi(RBI(op1), RBI(op2)));
2297 case N_FLONUM:
2298 return (cmp_flonum(bi_getd(RBI(op1)), RFF(op2)));
2299 case N_FIXRATIO:
2300 return (cmp_bi_fr(RBI(op1), RFRN(op2), RFRD(op2)));
2301 case N_BIGRATIO:
2302 return (cmp_bi_br(RBI(op1), RBR(op2)));
2303 }
2304 break;
2305 case N_FLONUM:
2306 switch (RTYPE(op2)) {
2307 case N_FIXNUM:
2308 return (cmp_flonum(RFF(op1), (double)RFI(op2)));
2309 case N_BIGNUM:
2310 return (cmp_flonum(RFF(op1), bi_getd(RBI(op2))));
2311 case N_FLONUM:
2312 return (cmp_flonum(RFF(op1), RFF(op2)));
2313 case N_FIXRATIO:
2314 return (cmp_flonum(RFF(op1),
2315 (double)RFRN(op2) / (double)RFRD(op2)));
2316 case N_BIGRATIO:
2317 return (cmp_flonum(RFF(op1), br_getd(RBR(op2))));
2318 }
2319 break;
2320 case N_FIXRATIO:
2321 switch (RTYPE(op2)) {
2322 case N_FIXNUM:
2323 return (cmp_fr_fi(RFRN(op1), RFRD(op1), RFI(op2)));
2324 case N_BIGNUM:
2325 return (cmp_fr_bi(RFRN(op1), RFRD(op1), RBI(op2)));
2326 case N_FLONUM:
2327 return (cmp_flonum((double)RFRN(op1) / (double)RFRD(op1),
2328 RFF(op2)));
2329 case N_FIXRATIO:
2330 return (cmp_fr_fr(RFRN(op1), RFRD(op1),
2331 RFRN(op2), RFRD(op2)));
2332 case N_BIGRATIO:
2333 return (cmp_fr_br(RFRN(op1), RFRD(op1), RBR(op2)));
2334 }
2335 break;
2336 case N_BIGRATIO:
2337 switch (RTYPE(op2)) {
2338 case N_FIXNUM:
2339 return (cmp_br_fi(RBR(op1), RFI(op2)));
2340 case N_BIGNUM:
2341 return (cmp_br_bi(RBR(op1), RBI(op2)));
2342 case N_FLONUM:
2343 return (cmp_flonum(br_getd(RBR(op1)), RFF(op2)));
2344 case N_FIXRATIO:
2345 return (cmp_br_fr(RBR(op1), RFRN(op2), RFRD(op2)));
2346 case N_BIGRATIO:
2347 return (cmp_br_br(RBR(op1), RBR(op2)));
2348 }
2349 }
2350
2351 return (0);
2352 }
2353
2354 static int
cmp_real_object(n_real * op1,LispObj * op2)2355 cmp_real_object(n_real *op1, LispObj *op2)
2356 {
2357 switch (OBJECT_TYPE(op2)) {
2358 case LispFixnum_t:
2359 switch (RTYPE(op1)) {
2360 case N_FIXNUM:
2361 return (cmp_fi_fi(RFI(op1), OFI(op2)));
2362 case N_BIGNUM:
2363 return (cmp_bi_fi(RBI(op1), OFI(op2)));
2364 case N_FLONUM:
2365 return (cmp_flonum(RFF(op1), (double)OFI(op2)));
2366 case N_FIXRATIO:
2367 return (cmp_fr_fi(RFRD(op1), RFRN(op1), OFI(op2)));
2368 case N_BIGRATIO:
2369 return (cmp_br_fi(RBR(op1), OFI(op2)));
2370 }
2371 break;
2372 case LispInteger_t:
2373 switch (RTYPE(op1)) {
2374 case N_FIXNUM:
2375 return (cmp_fi_fi(RFI(op1), OII(op2)));
2376 case N_BIGNUM:
2377 return (cmp_bi_fi(RBI(op1), OII(op2)));
2378 case N_FLONUM:
2379 return (cmp_flonum(RFF(op1), (double)OII(op2)));
2380 case N_FIXRATIO:
2381 return (cmp_fr_fi(RFRD(op1), RFRN(op1), OII(op2)));
2382 case N_BIGRATIO:
2383 return (cmp_br_fi(RBR(op1), OII(op2)));
2384 }
2385 break;
2386 case LispBignum_t:
2387 switch (RTYPE(op1)) {
2388 case N_FIXNUM:
2389 return (cmp_fi_bi(RFI(op1), OBI(op2)));
2390 case N_BIGNUM:
2391 return (cmp_bi_bi(RBI(op1), OBI(op2)));
2392 case N_FLONUM:
2393 return (cmp_flonum(RFF(op1), bi_getd(OBI(op2))));
2394 case N_FIXRATIO:
2395 return (cmp_fr_bi(RFRD(op1), RFRN(op1), OBI(op2)));
2396 case N_BIGRATIO:
2397 return (cmp_br_bi(RBR(op1), OBI(op2)));
2398 }
2399 break;
2400 case LispDFloat_t:
2401 switch (RTYPE(op1)) {
2402 case N_FIXNUM:
2403 return (cmp_flonum((double)RFI(op1), ODF(op2)));
2404 case N_BIGNUM:
2405 return (cmp_flonum(bi_getd(RBI(op1)), ODF(op2)));
2406 case N_FLONUM:
2407 return (cmp_flonum(RFF(op1), ODF(op2)));
2408 case N_FIXRATIO:
2409 return (cmp_flonum((double)RFRN(op1) / (double)RFRD(op1),
2410 ODF(op2)));
2411 case N_BIGRATIO:
2412 return (cmp_flonum(br_getd(RBR(op1)), ODF(op2)));
2413 }
2414 break;
2415 case LispRatio_t:
2416 switch (RTYPE(op1)) {
2417 case N_FIXNUM:
2418 return (cmp_fi_fr(RFI(op1), OFRN(op2), OFRD(op2)));
2419 case N_BIGNUM:
2420 return (cmp_bi_fr(RBI(op1), OFRN(op2), OFRD(op2)));
2421 case N_FLONUM:
2422 return (cmp_flonum(RFF(op1),
2423 (double)OFRN(op2) / (double)OFRD(op2)));
2424 case N_FIXRATIO:
2425 return (cmp_fr_fr(RFRN(op1), RFRD(op1),
2426 OFRN(op2), OFRD(op2)));
2427 case N_BIGRATIO:
2428 return (cmp_br_fr(RBR(op1), OFRN(op2), OFRD(op2)));
2429 }
2430 break;
2431 case LispBigratio_t:
2432 switch (RTYPE(op1)) {
2433 case N_FIXNUM:
2434 return (cmp_fi_br(RFI(op1), OBR(op2)));
2435 case N_BIGNUM:
2436 return (cmp_bi_br(RBI(op1), OBR(op2)));
2437 case N_FLONUM:
2438 return (cmp_flonum(RFF(op1), br_getd(OBR(op2))));
2439 case N_FIXRATIO:
2440 return (cmp_fr_br(RFRN(op1), RFRD(op1), OBR(op2)));
2441 case N_BIGRATIO:
2442 return (cmp_br_br(RBR(op1), OBR(op2)));
2443 }
2444 break;
2445 default:
2446 fatal_object_error(op2, NOT_A_REAL_NUMBER);
2447 break;
2448 }
2449
2450 return (0);
2451 }
2452
2453 #if 0 /* not used */
2454 static int
2455 cmp_number_object(n_number *op1, LispObj *op2)
2456 {
2457 if (op1->complex) {
2458 if (OBJECT_TYPE(op2) == LispComplex_t) {
2459 if (cmp_real_object(NREAL(op1), OCXR(op2)) == 0)
2460 return (cmp_real_object(NIMAG(op1), OCXI(op2)));
2461 return (1);
2462 }
2463 else if (cmp_real_real(NIMAG(op1), &zero) == 0)
2464 return (cmp_real_object(NREAL(op1), op2));
2465 else
2466 return (1);
2467 }
2468 else {
2469 switch (OBJECT_TYPE(op2)) {
2470 case LispFixnum_t:
2471 switch (NRTYPE(op1)) {
2472 case N_FIXNUM:
2473 return (cmp_fi_fi(NRFI(op1), OFI(op2)));
2474 case N_BIGNUM:
2475 return (cmp_bi_fi(NRBI(op1), OFI(op2)));
2476 case N_FLONUM:
2477 return (cmp_flonum(NRFF(op1), (double)OFI(op2)));
2478 case N_FIXRATIO:
2479 return (cmp_fr_fi(NRFRD(op1), NRFRN(op1), OFI(op2)));
2480 case N_BIGRATIO:
2481 return (cmp_br_fi(NRBR(op1), OFI(op2)));
2482 }
2483 break;
2484 case LispInteger_t:
2485 switch (NRTYPE(op1)) {
2486 case N_FIXNUM:
2487 return (cmp_fi_fi(NRFI(op1), OII(op2)));
2488 case N_BIGNUM:
2489 return (cmp_bi_fi(NRBI(op1), OII(op2)));
2490 case N_FLONUM:
2491 return (cmp_flonum(NRFF(op1), (double)OII(op2)));
2492 case N_FIXRATIO:
2493 return (cmp_fr_fi(NRFRD(op1), NRFRN(op1), OII(op2)));
2494 case N_BIGRATIO:
2495 return (cmp_br_fi(NRBR(op1), OII(op2)));
2496 }
2497 break;
2498 case LispBignum_t:
2499 switch (NRTYPE(op1)) {
2500 case N_FIXNUM:
2501 return (cmp_fi_bi(NRFI(op1), OBI(op2)));
2502 case N_BIGNUM:
2503 return (cmp_bi_bi(NRBI(op1), OBI(op2)));
2504 case N_FLONUM:
2505 return (cmp_flonum(NRFF(op1), bi_getd(OBI(op2))));
2506 case N_FIXRATIO:
2507 return (cmp_fr_bi(NRFRD(op1), NRFRN(op1), OBI(op2)));
2508 case N_BIGRATIO:
2509 return (cmp_br_bi(NRBR(op1), OBI(op2)));
2510 }
2511 break;
2512 case LispDFloat_t:
2513 switch (NRTYPE(op1)) {
2514 case N_FIXNUM:
2515 return (cmp_flonum((double)NRFI(op1), ODF(op2)));
2516 case N_BIGNUM:
2517 return (cmp_flonum(bi_getd(NRBI(op1)), ODF(op2)));
2518 case N_FLONUM:
2519 return (cmp_flonum(NRFF(op1), ODF(op2)));
2520 case N_FIXRATIO:
2521 return (cmp_flonum((double)NRFRN(op1) /
2522 (double)NRFRD(op1),
2523 ODF(op2)));
2524 case N_BIGRATIO:
2525 return (cmp_flonum(br_getd(NRBR(op1)), ODF(op2)));
2526 }
2527 break;
2528 case LispRatio_t:
2529 switch (NRTYPE(op1)) {
2530 case N_FIXNUM:
2531 return (cmp_fi_fr(NRFI(op1), OFRN(op2), OFRD(op2)));
2532 case N_BIGNUM:
2533 return (cmp_bi_fr(NRBI(op1), OFRN(op2), OFRD(op2)));
2534 case N_FLONUM:
2535 return (cmp_flonum(NRFF(op1),
2536 (double)OFRN(op2) / (double)OFRD(op2)));
2537 case N_FIXRATIO:
2538 return (cmp_fr_fr(NRFRN(op1), NRFRD(op1),
2539 OFRN(op2), OFRD(op2)));
2540 case N_BIGRATIO:
2541 return (cmp_br_fr(NRBR(op1), OFRN(op2), OFRD(op2)));
2542 }
2543 break;
2544 case LispBigratio_t:
2545 switch (NRTYPE(op1)) {
2546 case N_FIXNUM:
2547 return (cmp_fi_br(NRFI(op1), OBR(op2)));
2548 case N_BIGNUM:
2549 return (cmp_bi_br(NRBI(op1), OBR(op2)));
2550 case N_FLONUM:
2551 return (cmp_flonum(NRFF(op1), br_getd(OBR(op2))));
2552 case N_FIXRATIO:
2553 return (cmp_fr_br(NRFRN(op1), NRFRD(op1), OBR(op2)));
2554 case N_BIGRATIO:
2555 return (cmp_br_br(NRBR(op1), OBR(op2)));
2556 }
2557 break;
2558 case LispComplex_t:
2559 if (cmp_real_object(&zero, OCXI(op2)) == 0)
2560 return (cmp_real_object(NREAL(op1), OCXR(op2)));
2561 return (1);
2562 default:
2563 fatal_object_error(op2, NOT_A_NUMBER);
2564 break;
2565 }
2566 }
2567
2568 return (0);
2569 }
2570 #endif
2571
2572 static int
cmp_object_object(LispObj * op1,LispObj * op2,int real)2573 cmp_object_object(LispObj *op1, LispObj *op2, int real)
2574 {
2575 if (OBJECT_TYPE(op1) == LispComplex_t) {
2576 if (real)
2577 fatal_object_error(op1, NOT_A_REAL_NUMBER);
2578 if (OBJECT_TYPE(op2) == LispComplex_t)
2579 return (cmp_cx_cx(op1, op2));
2580 else if (cmp_real_object(&zero, OCXI(op1)) == 0)
2581 return (cmp_object_object(OCXR(op1), op2, real));
2582 return (1);
2583 }
2584 else if (OBJECT_TYPE(op2) == LispComplex_t) {
2585 if (real)
2586 fatal_object_error(op1, NOT_A_REAL_NUMBER);
2587 if (cmp_real_object(&zero, OCXI(op2)) == 0)
2588 return (cmp_object_object(op1, OCXR(op2), real));
2589 return (1);
2590 }
2591 else {
2592 switch (OBJECT_TYPE(op1)) {
2593 case LispFixnum_t:
2594 switch (OBJECT_TYPE(op2)) {
2595 case LispFixnum_t:
2596 return (cmp_fi_fi(OFI(op1), OFI(op2)));
2597 case LispInteger_t:
2598 return (cmp_fi_fi(OFI(op1), OII(op2)));
2599 case LispBignum_t:
2600 return (cmp_fi_bi(OFI(op1), OBI(op2)));
2601 case LispDFloat_t:
2602 return (cmp_flonum((double)OFI(op1), ODF(op2)));
2603 case LispRatio_t:
2604 return (cmp_fi_fr(OFI(op1),
2605 OFRN(op2), OFRD(op2)));
2606 case LispBigratio_t:
2607 return (cmp_fi_br(OFI(op1), OBR(op2)));
2608 default:
2609 break;
2610 }
2611 break;
2612 case LispInteger_t:
2613 switch (OBJECT_TYPE(op2)) {
2614 case LispFixnum_t:
2615 return (cmp_fi_fi(OII(op1), OFI(op2)));
2616 case LispInteger_t:
2617 return (cmp_fi_fi(OII(op1), OII(op2)));
2618 case LispBignum_t:
2619 return (cmp_fi_bi(OII(op1), OBI(op2)));
2620 case LispDFloat_t:
2621 return (cmp_flonum((double)OII(op1), ODF(op2)));
2622 case LispRatio_t:
2623 return (cmp_fi_fr(OII(op1),
2624 OFRN(op2), OFRD(op2)));
2625 case LispBigratio_t:
2626 return (cmp_fi_br(OII(op1), OBR(op2)));
2627 default:
2628 break;
2629 }
2630 break;
2631 case LispBignum_t:
2632 switch (OBJECT_TYPE(op2)) {
2633 case LispFixnum_t:
2634 return (cmp_bi_fi(OBI(op1), OFI(op2)));
2635 case LispInteger_t:
2636 return (cmp_bi_fi(OBI(op1), OII(op2)));
2637 case LispBignum_t:
2638 return (cmp_bi_bi(OBI(op1), OBI(op2)));
2639 case LispDFloat_t:
2640 return (cmp_flonum(bi_getd(OBI(op1)), ODF(op2)));
2641 case LispRatio_t:
2642 return (cmp_bi_fr(OBI(op1),
2643 OFRN(op2), OFRD(op2)));
2644 case LispBigratio_t:
2645 return (cmp_bi_br(OBI(op1), OBR(op2)));
2646 default:
2647 break;
2648 }
2649 break;
2650 case LispDFloat_t:
2651 switch (OBJECT_TYPE(op2)) {
2652 case LispFixnum_t:
2653 return (cmp_flonum(ODF(op1), (double)OFI(op2)));
2654 case LispInteger_t:
2655 return (cmp_flonum(ODF(op1), (double)OII(op2)));
2656 case LispBignum_t:
2657 return (cmp_flonum(ODF(op1), bi_getd(OBI(op2))));
2658 case LispDFloat_t:
2659 return (cmp_flonum(ODF(op1), ODF(op2)));
2660 break;
2661 case LispRatio_t:
2662 return (cmp_flonum(ODF(op1),
2663 (double)OFRN(op2) /
2664 (double)OFRD(op2)));
2665 case LispBigratio_t:
2666 return (cmp_flonum(ODF(op1), br_getd(OBR(op2))));
2667 default:
2668 break;
2669 }
2670 break;
2671 case LispRatio_t:
2672 switch (OBJECT_TYPE(op2)) {
2673 case LispFixnum_t:
2674 return (cmp_fr_fi(OFRN(op1), OFRD(op1), OFI(op2)));
2675 case LispInteger_t:
2676 return (cmp_fr_fi(OFRN(op1), OFRD(op1), OII(op2)));
2677 case LispBignum_t:
2678 return (cmp_fr_bi(OFRN(op1), OFRD(op1), OBI(op2)));
2679 case LispDFloat_t:
2680 return (cmp_flonum((double)OFRN(op1) /
2681 (double)OFRD(op1),
2682 ODF(op2)));
2683 case LispRatio_t:
2684 return (cmp_fr_fr(OFRN(op1), OFRD(op1),
2685 OFRN(op2), OFRD(op2)));
2686 case LispBigratio_t:
2687 return (cmp_fr_br(OFRN(op1), OFRD(op1), OBR(op2)));
2688 default:
2689 break;
2690 }
2691 break;
2692 case LispBigratio_t:
2693 switch (OBJECT_TYPE(op2)) {
2694 case LispFixnum_t:
2695 return (cmp_br_fi(OBR(op1), OFI(op2)));
2696 case LispInteger_t:
2697 return (cmp_br_fi(OBR(op1), OII(op2)));
2698 case LispBignum_t:
2699 return (cmp_br_bi(OBR(op1), OBI(op2)));
2700 case LispDFloat_t:
2701 return (cmp_flonum(br_getd(OBR(op1)), ODF(op2)));
2702 case LispRatio_t:
2703 return (cmp_br_fr(OBR(op1), OFRN(op2), OFRD(op2)));
2704 case LispBigratio_t:
2705 return (cmp_br_br(OBR(op1), OBR(op2)));
2706 default:
2707 break;
2708 }
2709 break;
2710 default:
2711 fatal_object_error(op1, NOT_A_NUMBER);
2712 break;
2713 }
2714 }
2715
2716 fatal_object_error(op2, NOT_A_NUMBER);
2717 return (0);
2718 }
2719
2720
2721 /************************************************************************
2722 * FIXNUM
2723 ************************************************************************/
2724 /*
2725 * check if op1 + op2 will overflow
2726 */
2727 static INLINE int
fi_fi_add_overflow(long op1,long op2)2728 fi_fi_add_overflow(long op1, long op2)
2729 {
2730 long op = op1 + op2;
2731
2732 return (op2 >= 0 ? op < op1 : op > op1);
2733 }
2734
2735 /*
2736 * check if op1 - op2 will overflow
2737 */
2738 static INLINE int
fi_fi_sub_overflow(long op1,long op2)2739 fi_fi_sub_overflow(long op1, long op2)
2740 {
2741 long op = op1 - op2;
2742
2743 return (op2 >= 0 ? op > op1 : op < op1);
2744 }
2745
2746 /*
2747 * check if op1 * op2 will overflow
2748 */
2749 static INLINE int
fi_fi_mul_overflow(long op1,long op2)2750 fi_fi_mul_overflow(long op1, long op2)
2751 {
2752 if (op1 == 0 || op1 == 1 || op2 == 0 || op2 == 1)
2753 return (0);
2754 if (op1 == MINSLONG || op2 == MINSLONG)
2755 return (1);
2756 if (op1 < 0)
2757 op1 = -op1;
2758 if (op2 < 0)
2759 op2 = -op2;
2760 return (op1 > MAXSLONG / op2);
2761 }
2762
2763
2764 /************************************************************************
2765 * BIGNUM
2766 ************************************************************************/
2767 static void
rbi_canonicalize(n_real * real)2768 rbi_canonicalize(n_real *real)
2769 {
2770 if (mpi_fiti(RBI(real))) {
2771 long fi = mpi_geti(RBI(real));
2772
2773 RTYPE(real) = N_FIXNUM;
2774 mpi_clear(RBI(real));
2775 XFREE(RBI(real));
2776 RFI(real) = fi;
2777 }
2778 }
2779
2780
2781 /************************************************************************
2782 * RATIO
2783 ************************************************************************/
2784 static void
rfr_canonicalize(n_real * real)2785 rfr_canonicalize(n_real *real)
2786 {
2787 long num, numerator, den, denominator, rest;
2788
2789 num = numerator = RFRN(real);
2790 den = denominator = RFRD(real);
2791 if (denominator == 0)
2792 fatal_error(DIVIDE_BY_ZERO);
2793
2794 if (num == MINSLONG || den == MINSLONG) {
2795 mpr *bigratio = XALLOC(mpr);
2796
2797 mpr_init(bigratio);
2798 mpr_seti(bigratio, num, den);
2799 RTYPE(real) = N_BIGRATIO;
2800 RBR(real) = bigratio;
2801 rbr_canonicalize(real);
2802 return;
2803 }
2804
2805 if (num < 0)
2806 num = -num;
2807 else if (num == 0) {
2808 RFI(real) = 0;
2809 RTYPE(real) = N_FIXNUM;
2810 return;
2811 }
2812 for (;;) {
2813 if ((rest = den % num) == 0)
2814 break;
2815 den = num;
2816 num = rest;
2817 }
2818 if (den != 1) {
2819 denominator /= num;
2820 numerator /= num;
2821 }
2822 if (denominator < 0) {
2823 numerator = -numerator;
2824 denominator = -denominator;
2825 }
2826 if (denominator == 1) {
2827 RTYPE(real) = N_FIXNUM;
2828 RFI(real) = numerator;
2829 }
2830 else {
2831 RFRN(real) = numerator;
2832 RFRD(real) = denominator;
2833 }
2834 }
2835
2836 static void
rbr_canonicalize(n_real * real)2837 rbr_canonicalize(n_real *real)
2838 {
2839 int fitnum, fitden;
2840 long numerator, denominator;
2841
2842 mpr_canonicalize(RBR(real));
2843 fitnum = mpi_fiti(RBRN(real));
2844 fitden = mpi_fiti(RBRD(real));
2845 if (fitnum && fitden) {
2846 numerator = mpi_geti(RBRN(real));
2847 denominator = mpi_geti(RBRD(real));
2848 mpr_clear(RBR(real));
2849 XFREE(RBR(real));
2850 if (numerator == 0) {
2851 RFI(real) = 0;
2852 RTYPE(real) = N_FIXNUM;
2853 }
2854 else if (denominator == 1) {
2855 RTYPE(real) = N_FIXNUM;
2856 RFI(real) = numerator;
2857 }
2858 else {
2859 RTYPE(real) = N_FIXRATIO;
2860 RFRN(real) = numerator;
2861 RFRD(real) = denominator;
2862 }
2863 }
2864 else if (fitden) {
2865 denominator = mpi_geti(RBRD(real));
2866 if (denominator == 1) {
2867 mpi *bigi = XALLOC(mpi);
2868
2869 mpi_init(bigi);
2870 mpi_set(bigi, RBRN(real));
2871 mpr_clear(RBR(real));
2872 XFREE(RBR(real));
2873 RTYPE(real) = N_BIGNUM;
2874 RBI(real) = bigi;
2875 }
2876 else if (denominator == 0)
2877 fatal_error(DIVIDE_BY_ZERO);
2878 }
2879 }
2880
2881
2882 /************************************************************************
2883 * COMPLEX
2884 ************************************************************************/
2885 static void
ncx_canonicalize(n_number * num)2886 ncx_canonicalize(n_number *num)
2887 {
2888 if (NITYPE(num) == N_FIXNUM && NIFI(num) == 0)
2889 num->complex = 0;
2890 }
2891
2892
2893 /************************************************************************
2894 * DIVIDE
2895 ************************************************************************/
2896 #define NDIVIDE_NOP 0
2897 #define NDIVIDE_ADD 1
2898 #define NDIVIDE_SUB 2
2899 static void
ndivide_fi_fi(n_number * num,long div,int fun,int flo)2900 ndivide_fi_fi(n_number *num, long div, int fun, int flo)
2901 {
2902 long quo, rem;
2903
2904 if (NRFI(num) == MINSLONG || div == MINSLONG) {
2905 LispObj integer;
2906 mpi *bignum = XALLOC(mpi);
2907
2908 mpi_init(bignum);
2909 mpi_seti(bignum, NRFI(num));
2910 NRBI(num) = bignum;
2911 NRTYPE(num) = N_BIGNUM;
2912 integer.type = LispInteger_t;
2913 integer.data.integer = div;
2914 ndivide_xi_xi(num, &integer, fun, flo);
2915 return;
2916 }
2917 else {
2918 quo = NRFI(num) / div;
2919 rem = NRFI(num) % div;
2920 }
2921
2922 switch (fun) {
2923 case NDIVIDE_CEIL:
2924 if ((rem < 0 && div < 0) || (rem > 0 && div > 0)) {
2925 ++quo;
2926 rem -= div;
2927 }
2928 break;
2929 case NDIVIDE_FLOOR:
2930 if ((rem < 0 && div > 0) || (rem > 0 && div < 0)) {
2931 --quo;
2932 rem += div;
2933 }
2934 break;
2935 case NDIVIDE_ROUND:
2936 if (div > 0) {
2937 if (rem > 0) {
2938 if (rem >= (div + 1) / 2) {
2939 ++quo;
2940 rem -= div;
2941 }
2942 }
2943 else {
2944 if (rem <= (-div - 1) / 2) {
2945 --quo;
2946 rem += div;
2947 }
2948 }
2949 }
2950 else {
2951 if (rem > 0) {
2952 if (rem >= (-div + 1) / 2) {
2953 --quo;
2954 rem += div;
2955 }
2956 }
2957 else {
2958 if (rem <= (div - 1) / 2) {
2959 ++quo;
2960 rem -= div;
2961 }
2962 }
2963 }
2964 break;
2965 }
2966
2967 NITYPE(num) = N_FIXNUM;
2968 NIFI(num) = rem;
2969 if (flo) {
2970 NRTYPE(num) = N_FLONUM;
2971 NRFF(num) = (double)quo;
2972 }
2973 else
2974 NRFI(num) = quo;
2975 }
2976
2977 static void
ndivide_xi_xi(n_number * num,LispObj * div,int fun,int flo)2978 ndivide_xi_xi(n_number *num, LispObj *div, int fun, int flo)
2979 {
2980 LispType type = OBJECT_TYPE(div);
2981 int state = NDIVIDE_NOP, dsign, rsign;
2982 mpi *quo, *rem;
2983
2984 quo = XALLOC(mpi);
2985 mpi_init(quo);
2986 if (NRTYPE(num) == N_FIXNUM)
2987 mpi_seti(quo, NRFI(num));
2988 else
2989 mpi_set(quo, NRBI(num));
2990
2991 rem = XALLOC(mpi);
2992 mpi_init(rem);
2993
2994 switch (type) {
2995 case LispFixnum_t:
2996 mpi_seti(rem, OFI(div));
2997 break;
2998 case LispInteger_t:
2999 mpi_seti(rem, OII(div));
3000 break;
3001 default:
3002 mpi_set(rem, OBI(div));
3003 }
3004
3005 dsign = mpi_sgn(rem);
3006
3007 mpi_divqr(quo, rem, quo, rem);
3008 rsign = mpi_sgn(rem);
3009
3010 switch (fun) {
3011 case NDIVIDE_CEIL:
3012 if ((rsign < 0 && dsign < 0) || (rsign > 0 && dsign > 0))
3013 state = NDIVIDE_ADD;
3014 break;
3015 case NDIVIDE_FLOOR:
3016 if ((rsign < 0 && dsign > 0) || (rsign > 0 && dsign < 0))
3017 state = NDIVIDE_SUB;
3018 break;
3019 case NDIVIDE_ROUND: {
3020 mpi test;
3021
3022 mpi_init(&test);
3023 switch (type) {
3024 case LispFixnum_t:
3025 mpi_seti(&test, OFI(div));
3026 break;
3027 case LispInteger_t:
3028 mpi_seti(&test, OII(div));
3029 break;
3030 default:
3031 mpi_set(&test, OBI(div));
3032 }
3033 if (dsign > 0) {
3034 if (rsign > 0) {
3035 mpi_addi(&test, &test, 1);
3036 mpi_divi(&test, &test, 2);
3037 if (mpi_cmp(rem, &test) >= 0)
3038 state = NDIVIDE_ADD;
3039 }
3040 else {
3041 mpi_neg(&test, &test);
3042 mpi_subi(&test, &test, 1);
3043 mpi_divi(&test, &test, 2);
3044 if (mpi_cmp(rem, &test) <= 0)
3045 state = NDIVIDE_SUB;
3046 }
3047 }
3048 else {
3049 if (rsign > 0) {
3050 mpi_neg(&test, &test);
3051 mpi_addi(&test, &test, 1);
3052 mpi_divi(&test, &test, 2);
3053 if (mpi_cmp(rem, &test) >= 0)
3054 state = NDIVIDE_SUB;
3055 }
3056 else {
3057 mpi_subi(&test, &test, 1);
3058 mpi_divi(&test, &test, 2);
3059 if (mpi_cmp(rem, &test) <= 0)
3060 state = NDIVIDE_ADD;
3061 }
3062 }
3063 mpi_clear(&test);
3064 } break;
3065 }
3066
3067 if (state == NDIVIDE_ADD) {
3068 mpi_addi(quo, quo, 1);
3069 switch (type) {
3070 case LispFixnum_t:
3071 mpi_subi(rem, rem, OFI(div));
3072 break;
3073 case LispInteger_t:
3074 mpi_subi(rem, rem, OII(div));
3075 break;
3076 default:
3077 mpi_sub(rem, rem, OBI(div));
3078 }
3079 }
3080 else if (state == NDIVIDE_SUB) {
3081 mpi_subi(quo, quo, 1);
3082 switch (type) {
3083 case LispFixnum_t:
3084 mpi_addi(rem, rem, OFI(div));
3085 break;
3086 case LispInteger_t:
3087 mpi_addi(rem, rem, OII(div));
3088 break;
3089 default:
3090 mpi_add(rem, rem, OBI(div));
3091 }
3092 }
3093
3094 if (mpi_fiti(rem)) {
3095 NITYPE(num) = N_FIXNUM;
3096 NIFI(num) = mpi_geti(rem);
3097 mpi_clear(rem);
3098 XFREE(rem);
3099 }
3100 else {
3101 NITYPE(num) = N_BIGNUM;
3102 NIBI(num) = rem;
3103 }
3104
3105 clear_real(NREAL(num));
3106
3107 if (flo) {
3108 double dval = bi_getd(quo);
3109
3110 mpi_clear(quo);
3111 XFREE(quo);
3112 NRTYPE(num) = N_FLONUM;
3113 NRFF(num) = dval;
3114 }
3115 else {
3116 NRTYPE(num) = N_BIGNUM;
3117 NRBI(num) = quo;
3118 rbi_canonicalize(NREAL(num));
3119 }
3120 }
3121
3122 static void
ndivide_flonum(n_number * number,double num,double div,int fun,int flo)3123 ndivide_flonum(n_number *number, double num, double div, int fun, int flo)
3124 {
3125 double quo, rem, modp, tmp;
3126
3127 modp = modf(num / div, &quo);
3128 rem = num - quo * div;
3129
3130 switch (fun) {
3131 case NDIVIDE_CEIL:
3132 if ((rem < 0.0 && div < 0.0) || (rem > 0.0 && div > 0.0)) {
3133 quo += 1.0;
3134 rem -= div;
3135 }
3136 break;
3137 case NDIVIDE_FLOOR:
3138 if ((rem < 0.0 && div > 0.0) || (rem > 0.0 && div < 0.0)) {
3139 quo -= 1.0;
3140 rem += div;
3141 }
3142 break;
3143 case NDIVIDE_ROUND:
3144 if (fabs(modp) != 0.5 || modf(quo * 0.5, &tmp) != 0.0) {
3145 if (div > 0.0) {
3146 if (rem > 0.0) {
3147 if (rem >= div * 0.5) {
3148 quo += 1.0;
3149 rem -= div;
3150 }
3151 }
3152 else {
3153 if (rem <= div * -0.5) {
3154 quo -= 1.0;
3155 rem += div;
3156 }
3157 }
3158 }
3159 else {
3160 if (rem > 0.0) {
3161 if (rem >= div * -0.5) {
3162 quo -= 1.0;
3163 rem += div;
3164 }
3165 }
3166 else {
3167 if (rem <= div * 0.5) {
3168 quo += 1.0;
3169 rem -= div;
3170 }
3171 }
3172 }
3173 }
3174 break;
3175 }
3176 if (!finite(quo) || !finite(rem))
3177 fatal_error(FLOATING_POINT_OVERFLOW);
3178
3179 NITYPE(number) = N_FLONUM;
3180 NIFF(number) = rem;
3181
3182 clear_real(NREAL(number));
3183
3184 if (flo) {
3185 NRTYPE(number) = N_FLONUM;
3186 NRFF(number) = quo;
3187 }
3188 else {
3189 if ((long)quo == quo) {
3190 NRTYPE(number) = N_FIXNUM;
3191 NRFI(number) = (long)quo;
3192 }
3193 else {
3194 mpi *bigi = XALLOC(mpi);
3195
3196 mpi_init(bigi);
3197 mpi_setd(bigi, quo);
3198 NRBI(number) = bigi;
3199 NRTYPE(number) = N_BIGNUM;
3200 }
3201 }
3202 }
3203
3204 static void
ndivide_xi_xr(n_number * num,LispObj * div,int fun,int flo)3205 ndivide_xi_xr(n_number *num, LispObj *div, int fun, int flo)
3206 {
3207 int state = NDIVIDE_NOP, dsign, rsign;
3208 mpi *quo;
3209 mpr *rem;
3210
3211 quo = XALLOC(mpi);
3212 mpi_init(quo);
3213 if (NRTYPE(num) == N_FIXNUM)
3214 mpi_seti(quo, NRFI(num));
3215 else
3216 mpi_set(quo, NRBI(num));
3217
3218 rem = XALLOC(mpr);
3219 mpr_init(rem);
3220
3221 if (XOBJECT_TYPE(div) == LispRatio_t)
3222 mpr_seti(rem, OFRN(div), OFRD(div));
3223 else
3224 mpr_set(rem, OBR(div));
3225 dsign = mpi_sgn(mpr_num(rem));
3226 mpi_mul(quo, quo, mpr_den(rem));
3227
3228 mpi_divqr(quo, mpr_num(rem), quo, mpr_num(rem));
3229 mpr_canonicalize(rem);
3230
3231 rsign = mpi_sgn(mpr_num(rem));
3232 if (mpr_fiti(rem)) {
3233 if (mpi_geti(mpr_den(rem)) == 1) {
3234 NITYPE(num) = N_FIXNUM;
3235 NIFI(num) = mpi_geti(mpr_num(rem));
3236 }
3237 else {
3238 NITYPE(num) = N_FIXRATIO;
3239 NIFRN(num) = mpi_geti(mpr_num(rem));
3240 NIFRD(num) = mpi_geti(mpr_den(rem));
3241 }
3242 mpr_clear(rem);
3243 XFREE(rem);
3244 }
3245 else {
3246 if (mpi_fiti(mpr_den(rem)) && mpi_geti(mpr_den(rem)) == 1) {
3247 NITYPE(num) = N_BIGNUM;
3248 NIBI(num) = mpr_num(rem);
3249 mpi_clear(mpr_den(rem));
3250 XFREE(rem);
3251 }
3252 else {
3253 NITYPE(num) = N_BIGRATIO;
3254 NIBR(num) = rem;
3255 }
3256 }
3257
3258 switch (fun) {
3259 case NDIVIDE_CEIL:
3260 if ((rsign < 0 && dsign < 0) || (rsign > 0 && dsign > 0))
3261 state = NDIVIDE_ADD;
3262 break;
3263 case NDIVIDE_FLOOR:
3264 if ((rsign < 0 && dsign > 0) || (rsign > 0 && dsign < 0))
3265 state = NDIVIDE_SUB;
3266 break;
3267 case NDIVIDE_ROUND: {
3268 n_real cmp;
3269
3270 set_real_object(&cmp, div);
3271 div_real_real(&cmp, &two);
3272 if (dsign > 0) {
3273 if (rsign > 0) {
3274 if (cmp_real_real(NIMAG(num), &cmp) >= 0)
3275 state = NDIVIDE_ADD;
3276 }
3277 else {
3278 neg_real(&cmp);
3279 if (cmp_real_real(NIMAG(num), &cmp) <= 0)
3280 state = NDIVIDE_SUB;
3281 }
3282 }
3283 else {
3284 if (rsign > 0) {
3285 neg_real(&cmp);
3286 if (cmp_real_real(NIMAG(num), &cmp) >= 0)
3287 state = NDIVIDE_SUB;
3288 }
3289 else {
3290 if (cmp_real_real(NIMAG(num), &cmp) <= 0)
3291 state = NDIVIDE_ADD;
3292 }
3293 }
3294 clear_real(&cmp);
3295 } break;
3296 }
3297
3298 if (state == NDIVIDE_ADD) {
3299 mpi_addi(quo, quo, 1);
3300 sub_real_object(NIMAG(num), div);
3301 }
3302 else if (state == NDIVIDE_SUB) {
3303 mpi_subi(quo, quo, 1);
3304 add_real_object(NIMAG(num), div);
3305 }
3306
3307 clear_real(NREAL(num));
3308
3309 if (flo) {
3310 double dval = bi_getd(quo);
3311
3312 mpi_clear(quo);
3313 XFREE(quo);
3314 NRTYPE(num) = N_FLONUM;
3315 NRFF(num) = dval;
3316 }
3317 else {
3318 NRBI(num) = quo;
3319 NRTYPE(num) = N_BIGNUM;
3320 rbi_canonicalize(NREAL(num));
3321 }
3322 }
3323
3324 static void
ndivide_xr_xi(n_number * num,LispObj * div,int fun,int flo)3325 ndivide_xr_xi(n_number *num, LispObj *div, int fun, int flo)
3326 {
3327 LispType type = OBJECT_TYPE(div);
3328 int state = NDIVIDE_NOP, dsign, rsign;
3329 mpi *quo;
3330 mpr *rem;
3331
3332 quo = XALLOC(mpi);
3333 mpi_init(quo);
3334 switch (type) {
3335 case LispFixnum_t:
3336 dsign = OFI(div) < 0 ? -1 : OFI(div) > 0 ? 1 : 0;
3337 mpi_seti(quo, OFI(div));
3338 break;
3339 case LispInteger_t:
3340 dsign = OII(div) < 0 ? -1 : OII(div) > 0 ? 1 : 0;
3341 mpi_seti(quo, OII(div));
3342 break;
3343 default:
3344 dsign = mpi_sgn(OBI(div));
3345 mpi_set(quo, OBI(div));
3346 break;
3347 }
3348
3349 rem = XALLOC(mpr);
3350 mpr_init(rem);
3351 if (NRTYPE(num) == N_FIXRATIO) {
3352 mpr_seti(rem, NRFRN(num), NRFRD(num));
3353 mpi_muli(quo, quo, NRFRD(num));
3354 }
3355 else {
3356 mpr_set(rem, NRBR(num));
3357 mpi_mul(quo, quo, NRBRD(num));
3358 }
3359 mpi_divqr(quo, mpr_num(rem), mpr_num(rem), quo);
3360 mpr_canonicalize(rem);
3361
3362 rsign = mpi_sgn(mpr_num(rem));
3363 if (mpr_fiti(rem)) {
3364 NITYPE(num) = N_FIXRATIO;
3365 NIFRN(num) = mpi_geti(mpr_num(rem));
3366 NIFRD(num) = mpi_geti(mpr_den(rem));
3367 mpr_clear(rem);
3368 XFREE(rem);
3369 }
3370 else {
3371 NITYPE(num) = N_BIGRATIO;
3372 NIBR(num) = rem;
3373 }
3374
3375 switch (fun) {
3376 case NDIVIDE_CEIL:
3377 if ((rsign < 0 && dsign < 0) || (rsign > 0 && dsign > 0))
3378 state = NDIVIDE_ADD;
3379 break;
3380 case NDIVIDE_FLOOR:
3381 if ((rsign < 0 && dsign > 0) || (rsign > 0 && dsign < 0))
3382 state = NDIVIDE_SUB;
3383 break;
3384 case NDIVIDE_ROUND: {
3385 n_real cmp;
3386
3387 set_real_object(&cmp, div);
3388 div_real_real(&cmp, &two);
3389 if (dsign > 0) {
3390 if (rsign > 0) {
3391 if (cmp_real_real(NIMAG(num), &cmp) >= 0)
3392 state = NDIVIDE_ADD;
3393 }
3394 else {
3395 neg_real(&cmp);
3396 if (cmp_real_real(NIMAG(num), &cmp) <= 0)
3397 state = NDIVIDE_SUB;
3398 }
3399 }
3400 else {
3401 if (rsign > 0) {
3402 neg_real(&cmp);
3403 if (cmp_real_real(NIMAG(num), &cmp) >= 0)
3404 state = NDIVIDE_SUB;
3405 }
3406 else {
3407 if (cmp_real_real(NIMAG(num), &cmp) <= 0)
3408 state = NDIVIDE_ADD;
3409 }
3410 }
3411 clear_real(&cmp);
3412 } break;
3413 }
3414
3415 if (state == NDIVIDE_ADD) {
3416 mpi_addi(quo, quo, 1);
3417 sub_real_object(NIMAG(num), div);
3418 }
3419 else if (state == NDIVIDE_SUB) {
3420 mpi_subi(quo, quo, 1);
3421 add_real_object(NIMAG(num), div);
3422 }
3423
3424 clear_real(NREAL(num));
3425
3426 if (flo) {
3427 double dval = bi_getd(quo);
3428
3429 mpi_clear(quo);
3430 XFREE(quo);
3431 NRTYPE(num) = N_FLONUM;
3432 NRFF(num) = dval;
3433 }
3434 else {
3435 NRBI(num) = quo;
3436 NRTYPE(num) = N_BIGNUM;
3437 rbi_canonicalize(NREAL(num));
3438 }
3439 }
3440
3441 static void
ndivide_xr_xr(n_number * num,LispObj * div,int fun,int flo)3442 ndivide_xr_xr(n_number *num, LispObj *div, int fun, int flo)
3443 {
3444 int state = NDIVIDE_NOP, dsign, rsign, modp;
3445 mpr *bigr;
3446 mpi *bigi;
3447
3448 bigr = XALLOC(mpr);
3449 mpr_init(bigr);
3450 if (NRTYPE(num) == N_FIXRATIO)
3451 mpr_seti(bigr, NRFRN(num), NRFRD(num));
3452 else
3453 mpr_set(bigr, NRBR(num));
3454
3455 NITYPE(num) = N_BIGRATIO;
3456 NIBR(num) = bigr;
3457
3458 if (OBJECT_TYPE(div) == LispRatio_t) {
3459 dsign = OFRN(div) < 0 ? -1 : OFRN(div) > 0 ? 1 : 0;
3460 mpi_muli(mpr_num(bigr), mpr_num(bigr), OFRD(div));
3461 mpi_muli(mpr_den(bigr), mpr_den(bigr), OFRN(div));
3462 }
3463 else {
3464 dsign = mpi_sgn(OBRN(div));
3465 mpr_div(bigr, bigr, OBR(div));
3466 }
3467 modp = mpi_fiti(mpr_den(bigr)) && mpi_geti(mpr_den(bigr)) == 2;
3468
3469 bigi = XALLOC(mpi);
3470 mpi_init(bigi);
3471 mpi_divqr(bigi, mpr_num(bigr), mpr_num(bigr), mpr_den(bigr));
3472
3473 if (OBJECT_TYPE(div) == LispRatio_t)
3474 mpi_seti(mpr_den(bigr), OFRD(div));
3475 else
3476 mpi_set(mpr_den(bigr), OBRD(div));
3477 if (NRTYPE(num) == N_FIXRATIO)
3478 mpi_muli(mpr_den(bigr), mpr_den(bigr), NRFRD(num));
3479 else
3480 mpi_mul(mpr_den(bigr), mpr_den(bigr), NRBRD(num));
3481
3482 clear_real(NREAL(num));
3483 NRTYPE(num) = N_BIGNUM;
3484 NRBI(num) = bigi;
3485
3486 rbr_canonicalize(NIMAG(num));
3487 rsign = cmp_real_real(NIMAG(num), &zero);
3488
3489 switch (fun) {
3490 case NDIVIDE_CEIL:
3491 if ((rsign < 0 && dsign < 0) || (rsign > 0 && dsign > 0))
3492 state = NDIVIDE_ADD;
3493 break;
3494 case NDIVIDE_FLOOR:
3495 if ((rsign < 0 && dsign > 0) || (rsign > 0 && dsign < 0))
3496 state = NDIVIDE_SUB;
3497 break;
3498 case NDIVIDE_ROUND:
3499 if (!modp || (bigi->digs[0] & 1) == 1) {
3500 n_real cmp;
3501
3502 set_real_object(&cmp, div);
3503 div_real_real(&cmp, &two);
3504 if (dsign > 0) {
3505 if (rsign > 0) {
3506 if (cmp_real_real(NIMAG(num), &cmp) >= 0)
3507 state = NDIVIDE_ADD;
3508 }
3509 else {
3510 neg_real(&cmp);
3511 if (cmp_real_real(NIMAG(num), &cmp) <= 0)
3512 state = NDIVIDE_SUB;
3513 }
3514 }
3515 else {
3516 if (rsign > 0) {
3517 neg_real(&cmp);
3518 if (cmp_real_real(NIMAG(num), &cmp) >= 0)
3519 state = NDIVIDE_SUB;
3520 }
3521 else {
3522 if (cmp_real_real(NIMAG(num), &cmp) <= 0)
3523 state = NDIVIDE_ADD;
3524 }
3525 }
3526 clear_real(&cmp);
3527 }
3528 break;
3529 }
3530
3531 if (state == NDIVIDE_ADD) {
3532 add_real_real(NREAL(num), &one);
3533 sub_real_object(NIMAG(num), div);
3534 }
3535 else if (state == NDIVIDE_SUB) {
3536 sub_real_real(NREAL(num), &one);
3537 add_real_object(NIMAG(num), div);
3538 }
3539
3540 if (NRTYPE(num) == N_BIGNUM) {
3541 if (flo) {
3542 double dval = bi_getd(bigi);
3543
3544 mpi_clear(bigi);
3545 XFREE(bigi);
3546 NRTYPE(num) = N_FLONUM;
3547 NRFF(num) = dval;
3548 }
3549 else
3550 rbi_canonicalize(NREAL(num));
3551 }
3552 else if (flo) {
3553 NRTYPE(num) = N_FLONUM;
3554 NRFF(num) = (double)NRFI(num);
3555 }
3556 }
3557
3558
3559 /************************************************************************
3560 * REAL COMPLEX
3561 ************************************************************************/
3562 static void
nadd_re_cx(n_number * num,LispObj * comp)3563 nadd_re_cx(n_number *num, LispObj *comp)
3564 {
3565 /*
3566 Ra+Rb Ib
3567 */
3568 /* Ra+Rb */
3569 add_real_object(NREAL(num), OCXR(comp));
3570
3571 /* Ib */
3572 set_real_object(NIMAG(num), OCXI(comp));
3573
3574 num->complex = 1;
3575
3576 ncx_canonicalize(num);
3577 }
3578
3579 static void
nsub_re_cx(n_number * num,LispObj * comp)3580 nsub_re_cx(n_number *num, LispObj *comp)
3581 {
3582 /*
3583 Ra-Rb -Ib
3584 */
3585 /* Ra-Rb */
3586 sub_real_object(NREAL(num), OCXR(comp));
3587
3588 /* -Ib */
3589 NITYPE(num) = N_FIXNUM;
3590 NIFI(num) = -1;
3591 mul_real_object(NIMAG(num), OCXI(comp));
3592
3593 num->complex = 1;
3594
3595 ncx_canonicalize(num);
3596 }
3597
3598 static void
nmul_re_cx(n_number * num,LispObj * comp)3599 nmul_re_cx(n_number *num, LispObj *comp)
3600 {
3601 /*
3602 Ra*Rb Ra*Ib
3603 */
3604 /* copy before change */
3605 set_real_real(NIMAG(num), NREAL(num));
3606
3607 /* Ra*Rb */
3608 mul_real_object(NREAL(num), OCXR(comp));
3609
3610 /* Ra*Ib */
3611 mul_real_object(NIMAG(num), OCXI(comp));
3612
3613 num->complex = 1;
3614
3615 ncx_canonicalize(num);
3616 }
3617
3618 static void
ndiv_re_cx(n_number * num,LispObj * comp)3619 ndiv_re_cx(n_number *num, LispObj *comp)
3620 {
3621 /*
3622 Ra*Rb -Ib*Ra
3623 ----------- -----------
3624 Rb*Rb+Ib*Ib Rb*Rb+Ib*Ib
3625 */
3626 n_real div, temp;
3627
3628 /* Rb*Rb */
3629 set_real_object(&div, OCXR(comp));
3630 mul_real_object(&div, OCXR(comp));
3631
3632 /* Ib*Ib */
3633 set_real_object(&temp, OCXI(comp));
3634 mul_real_object(&temp, OCXI(comp));
3635
3636 /* Rb*Rb+Ib*Ib */
3637 add_real_real(&div, &temp);
3638 clear_real(&temp);
3639
3640 /* -Ib*Ra */
3641 NITYPE(num) = N_FIXNUM;
3642 NIFI(num) = -1;
3643 mul_real_object(NIMAG(num), OCXI(comp));
3644 mul_real_real(NIMAG(num), NREAL(num));
3645
3646 /* Ra*Rb */
3647 mul_real_object(NREAL(num), OCXR(comp));
3648
3649 div_real_real(NREAL(num), &div);
3650 div_real_real(NIMAG(num), &div);
3651 clear_real(&div);
3652
3653 num->complex = 1;
3654
3655 ncx_canonicalize(num);
3656 }
3657
3658
3659 /************************************************************************
3660 * COMPLEX REAL
3661 ************************************************************************/
3662 static void
nadd_cx_re(n_number * num,LispObj * re)3663 nadd_cx_re(n_number *num, LispObj *re)
3664 {
3665 /*
3666 Ra+Rb Ia
3667 */
3668 add_real_object(NREAL(num), re);
3669
3670 ncx_canonicalize(num);
3671 }
3672
3673 static void
nsub_cx_re(n_number * num,LispObj * re)3674 nsub_cx_re(n_number *num, LispObj *re)
3675 {
3676 /*
3677 Ra-Rb Ia
3678 */
3679 sub_real_object(NREAL(num), re);
3680
3681 ncx_canonicalize(num);
3682 }
3683
3684 static void
nmul_cx_re(n_number * num,LispObj * re)3685 nmul_cx_re(n_number *num, LispObj *re)
3686 {
3687 /*
3688 Ra*Rb Ia*Rb
3689 */
3690 mul_real_object(NREAL(num), re);
3691 mul_real_object(NIMAG(num), re);
3692
3693 ncx_canonicalize(num);
3694 }
3695
3696 static void
ndiv_cx_re(n_number * num,LispObj * re)3697 ndiv_cx_re(n_number *num, LispObj *re)
3698 {
3699 /*
3700 Ra/Rb Ia/Rb
3701 */
3702 div_real_object(NREAL(num), re);
3703 div_real_object(NIMAG(num), re);
3704
3705 ncx_canonicalize(num);
3706 }
3707
3708
3709 /************************************************************************
3710 * COMPLEX COMPLEX
3711 ************************************************************************/
3712 static void
nadd_cx_cx(n_number * num,LispObj * comp)3713 nadd_cx_cx(n_number *num, LispObj *comp)
3714 {
3715 /*
3716 Ra+Rb Ia+Ib
3717 */
3718 add_real_object(NREAL(num), OCXR(comp));
3719 add_real_object(NIMAG(num), OCXI(comp));
3720
3721 ncx_canonicalize(num);
3722 }
3723
3724 static void
nsub_cx_cx(n_number * num,LispObj * comp)3725 nsub_cx_cx(n_number *num, LispObj *comp)
3726 {
3727 /*
3728 Ra-Rb Ia-Ib
3729 */
3730 sub_real_object(NREAL(num), OCXR(comp));
3731 sub_real_object(NIMAG(num), OCXI(comp));
3732
3733 ncx_canonicalize(num);
3734 }
3735
3736 static void
nmul_cx_cx(n_number * num,LispObj * comp)3737 nmul_cx_cx(n_number *num, LispObj *comp)
3738 {
3739 /*
3740 Ra*Rb-Ia*Ib Ra*Ib+Ia*Rb
3741 */
3742 n_real IaIb, RaIb;
3743
3744 set_real_real(&IaIb, NIMAG(num));
3745 mul_real_object(&IaIb, OCXI(comp));
3746
3747 set_real_real(&RaIb, NREAL(num));
3748 mul_real_object(&RaIb, OCXI(comp));
3749
3750 /* Ra*Rb-Ia*Ib */
3751 mul_real_object(NREAL(num), OCXR(comp));
3752 sub_real_real(NREAL(num), &IaIb);
3753 clear_real(&IaIb);
3754
3755 /* Ra*Ib+Ia*Rb */
3756 mul_real_object(NIMAG(num), OCXR(comp));
3757 add_real_real(NIMAG(num), &RaIb);
3758 clear_real(&RaIb);
3759
3760 ncx_canonicalize(num);
3761 }
3762
3763 static void
ndiv_cx_cx(n_number * num,LispObj * comp)3764 ndiv_cx_cx(n_number *num, LispObj *comp)
3765 {
3766 /*
3767 Ra*Rb+Ia*Ib Ia*Rb-Ib*Ra
3768 ----------- -----------
3769 Rb*Rb+Ib*Ib Rb*Rb+Ib*Ib
3770 */
3771 n_real temp1, temp2;
3772
3773 /* IaIb */
3774 set_real_real(&temp1, NIMAG(num));
3775 mul_real_object(&temp1, OCXI(comp));
3776
3777 /* IbRa */
3778 set_real_real(&temp2, NREAL(num));
3779 mul_real_object(&temp2, OCXI(comp));
3780
3781 /* Ra*Rb+Ia*Ib */
3782 mul_real_object(NREAL(num), OCXR(comp));
3783 add_real_real(NREAL(num), &temp1);
3784 clear_real(&temp1);
3785
3786 /* Ia*Rb-Ib*Ra */
3787 mul_real_object(NIMAG(num), OCXR(comp));
3788 sub_real_real(NIMAG(num), &temp2);
3789 clear_real(&temp2);
3790
3791
3792 /* Rb*Rb */
3793 set_real_object(&temp1, OCXR(comp));
3794 mul_real_object(&temp1, OCXR(comp));
3795
3796 /* Ib*Ib */
3797 set_real_object(&temp2, OCXI(comp));
3798 mul_real_object(&temp2, OCXI(comp));
3799
3800 /* Rb*Rb+Ib*Ib */
3801 add_real_real(&temp1, &temp2);
3802 clear_real(&temp2);
3803
3804 div_real_real(NREAL(num), &temp1);
3805 div_real_real(NIMAG(num), &temp1);
3806 clear_real(&temp1);
3807
3808 ncx_canonicalize(num);
3809 }
3810
3811 static int
cmp_cx_cx(LispObj * op1,LispObj * op2)3812 cmp_cx_cx(LispObj *op1, LispObj *op2)
3813 {
3814 int cmp;
3815
3816 cmp = cmp_object_object(OCXR(op1), OCXR(op2), 1);
3817 if (cmp == 0)
3818 cmp = cmp_object_object(OCXI(op1), OCXI(op2), 1);
3819
3820 return (cmp);
3821 }
3822
3823
3824 /************************************************************************
3825 * FLONUM FLONUM
3826 ************************************************************************/
3827 static void
radd_flonum(n_real * real,double op1,double op2)3828 radd_flonum(n_real *real, double op1, double op2)
3829 {
3830 double value = op1 + op2;
3831
3832 if (!finite(value))
3833 fatal_error(FLOATING_POINT_OVERFLOW);
3834 switch (RTYPE(real)) {
3835 case N_FIXNUM:
3836 case N_FIXRATIO:
3837 RTYPE(real) = N_FLONUM;
3838 break;
3839 case N_BIGNUM:
3840 RCLEAR_BI(real);
3841 RTYPE(real) = N_FLONUM;
3842 break;
3843 case N_BIGRATIO:
3844 RCLEAR_BR(real);
3845 RTYPE(real) = N_FLONUM;
3846 break;
3847 }
3848 RFF(real) = value;
3849 }
3850
3851 static void
rsub_flonum(n_real * real,double op1,double op2)3852 rsub_flonum(n_real *real, double op1, double op2)
3853 {
3854 double value = op1 - op2;
3855
3856 if (!finite(value))
3857 fatal_error(FLOATING_POINT_OVERFLOW);
3858 switch (RTYPE(real)) {
3859 case N_FIXNUM:
3860 case N_FIXRATIO:
3861 RTYPE(real) = N_FLONUM;
3862 break;
3863 case N_BIGNUM:
3864 RCLEAR_BI(real);
3865 RTYPE(real) = N_FLONUM;
3866 break;
3867 case N_BIGRATIO:
3868 RCLEAR_BR(real);
3869 RTYPE(real) = N_FLONUM;
3870 break;
3871 }
3872 RFF(real) = value;
3873 }
3874
3875 static void
rmul_flonum(n_real * real,double op1,double op2)3876 rmul_flonum(n_real *real, double op1, double op2)
3877 {
3878 double value = op1 * op2;
3879
3880 if (!finite(value))
3881 fatal_error(FLOATING_POINT_OVERFLOW);
3882 switch (RTYPE(real)) {
3883 case N_FIXNUM:
3884 case N_FIXRATIO:
3885 RTYPE(real) = N_FLONUM;
3886 break;
3887 case N_BIGNUM:
3888 RCLEAR_BI(real);
3889 RTYPE(real) = N_FLONUM;
3890 break;
3891 case N_BIGRATIO:
3892 RCLEAR_BR(real);
3893 RTYPE(real) = N_FLONUM;
3894 break;
3895 }
3896 RFF(real) = value;
3897 }
3898
3899 static void
rdiv_flonum(n_real * real,double op1,double op2)3900 rdiv_flonum(n_real *real, double op1, double op2)
3901 {
3902 double value;
3903
3904 if (op2 == 0.0)
3905 fatal_error(DIVIDE_BY_ZERO);
3906 value = op1 / op2;
3907 if (!finite(value))
3908 fatal_error(FLOATING_POINT_OVERFLOW);
3909 switch (RTYPE(real)) {
3910 case N_FIXNUM:
3911 case N_FIXRATIO:
3912 RTYPE(real) = N_FLONUM;
3913 break;
3914 case N_BIGNUM:
3915 RCLEAR_BI(real);
3916 RTYPE(real) = N_FLONUM;
3917 break;
3918 case N_BIGRATIO:
3919 RCLEAR_BR(real);
3920 RTYPE(real) = N_FLONUM;
3921 break;
3922 }
3923 RFF(real) = value;
3924 }
3925
3926 static int
cmp_flonum(double op1,double op2)3927 cmp_flonum(double op1, double op2)
3928 {
3929 double value = op1 - op2;
3930
3931 if (!finite(value))
3932 fatal_error(FLOATING_POINT_OVERFLOW);
3933
3934 return (value > 0.0 ? 1 : value < 0.0 ? -1 : 0);
3935 }
3936
3937
3938 /************************************************************************
3939 * FIXNUM FIXNUM
3940 ************************************************************************/
3941 static void
rop_fi_fi_bi(n_real * real,long fi,int op)3942 rop_fi_fi_bi(n_real *real, long fi, int op)
3943 {
3944 mpi *bigi = XALLOC(mpi);
3945
3946 mpi_init(bigi);
3947 mpi_seti(bigi, RFI(real));
3948 if (op == NOP_ADD)
3949 mpi_addi(bigi, bigi, fi);
3950 else if (op == NOP_SUB)
3951 mpi_subi(bigi, bigi, fi);
3952 else
3953 mpi_muli(bigi, bigi, fi);
3954 RBI(real) = bigi;
3955 RTYPE(real) = N_BIGNUM;
3956 }
3957
3958 static INLINE void
radd_fi_fi(n_real * real,long fi)3959 radd_fi_fi(n_real *real, long fi)
3960 {
3961 if (!fi_fi_add_overflow(RFI(real), fi))
3962 RFI(real) += fi;
3963 else
3964 rop_fi_fi_bi(real, fi, NOP_ADD);
3965 }
3966
3967 static INLINE void
rsub_fi_fi(n_real * real,long fi)3968 rsub_fi_fi(n_real *real, long fi)
3969 {
3970 if (!fi_fi_sub_overflow(RFI(real), fi))
3971 RFI(real) -= fi;
3972 else
3973 rop_fi_fi_bi(real, fi, NOP_SUB);
3974 }
3975
3976 static INLINE void
rmul_fi_fi(n_real * real,long fi)3977 rmul_fi_fi(n_real *real, long fi)
3978 {
3979 if (!fi_fi_mul_overflow(RFI(real), fi))
3980 RFI(real) *= fi;
3981 else
3982 rop_fi_fi_bi(real, fi, NOP_MUL);
3983 }
3984
3985 static INLINE void
rdiv_fi_fi(n_real * real,long fi)3986 rdiv_fi_fi(n_real *real, long fi)
3987 {
3988 RTYPE(real) = N_FIXRATIO;
3989 RFRN(real) = RFI(real);
3990 RFRD(real) = fi;
3991 rfr_canonicalize(real);
3992 }
3993
3994 static INLINE int
cmp_fi_fi(long op1,long op2)3995 cmp_fi_fi(long op1, long op2)
3996 {
3997 if (op1 > op2)
3998 return (1);
3999 else if (op1 < op2)
4000 return (-1);
4001
4002 return (0);
4003 }
4004
4005
4006 /************************************************************************
4007 * FIXNUM BIGNUM
4008 ************************************************************************/
4009 static void
rop_fi_bi_xi(n_real * real,mpi * bi,int nop)4010 rop_fi_bi_xi(n_real *real, mpi *bi, int nop)
4011 {
4012 mpi *bigi = XALLOC(mpi);
4013
4014 mpi_init(bigi);
4015 mpi_seti(bigi, RFI(real));
4016 if (nop == NOP_ADD)
4017 mpi_add(bigi, bigi, bi);
4018 else if (nop == NOP_SUB)
4019 mpi_sub(bigi, bigi, bi);
4020 else
4021 mpi_mul(bigi, bigi, bi);
4022
4023 if (mpi_fiti(bigi)) {
4024 RFI(real) = mpi_geti(bigi);
4025 mpi_clear(bigi);
4026 XFREE(bigi);
4027 }
4028 else {
4029 RBI(real) = bigi;
4030 RTYPE(real) = N_BIGNUM;
4031 }
4032 }
4033
4034 static INLINE void
radd_fi_bi(n_real * real,mpi * bi)4035 radd_fi_bi(n_real *real, mpi *bi)
4036 {
4037 rop_fi_bi_xi(real, bi, NOP_ADD);
4038 }
4039
4040 static INLINE void
rsub_fi_bi(n_real * real,mpi * bi)4041 rsub_fi_bi(n_real *real, mpi *bi)
4042 {
4043 rop_fi_bi_xi(real, bi, NOP_SUB);
4044 }
4045
4046 static INLINE void
rmul_fi_bi(n_real * real,mpi * bi)4047 rmul_fi_bi(n_real *real, mpi *bi)
4048 {
4049 rop_fi_bi_xi(real, bi, NOP_MUL);
4050 }
4051
4052 static void
rdiv_fi_bi(n_real * real,mpi * bi)4053 rdiv_fi_bi(n_real *real, mpi *bi)
4054 {
4055 mpr *bigr;
4056
4057 if (mpi_cmpi(bi, 0) == 0)
4058 fatal_error(DIVIDE_BY_ZERO);
4059
4060 bigr = XALLOC(mpr);
4061 mpr_init(bigr);
4062 mpi_seti(mpr_num(bigr), RFI(real));
4063 mpi_set(mpr_den(bigr), bi);
4064 RBR(real) = bigr;
4065 RTYPE(real) = N_BIGRATIO;
4066 rbr_canonicalize(real);
4067 }
4068
4069 static INLINE int
cmp_fi_bi(long fixnum,mpi * bignum)4070 cmp_fi_bi(long fixnum, mpi *bignum)
4071 {
4072 return (-mpi_cmpi(bignum, fixnum));
4073 }
4074
4075
4076 /************************************************************************
4077 * FIXNUM FIXRATIO
4078 ************************************************************************/
4079 static void
rop_fi_fr_as_xr(n_real * real,long num,long den,int nop)4080 rop_fi_fr_as_xr(n_real *real, long num, long den, int nop)
4081 {
4082 int fit;
4083 long value = 0, op = RFI(real);
4084
4085 fit = !fi_fi_mul_overflow(op, den);
4086 if (fit) {
4087 value = op * den;
4088 if (nop == NOP_ADD)
4089 fit = !fi_fi_add_overflow(value, num);
4090 else
4091 fit = !fi_fi_sub_overflow(value, num);
4092 }
4093 if (fit) {
4094 if (nop == NOP_ADD)
4095 RFRN(real) = value + num;
4096 else
4097 RFRN(real) = value - num;
4098 RFRD(real) = den;
4099 RTYPE(real) = N_FIXRATIO;
4100 rfr_canonicalize(real);
4101 }
4102 else {
4103 mpi iop;
4104 mpr *bigr = XALLOC(mpr);
4105
4106 mpi_init(&iop);
4107 mpi_seti(&iop, op);
4108 mpi_muli(&iop, &iop, den);
4109
4110 mpr_init(bigr);
4111 mpr_seti(bigr, num, den);
4112 if (nop == NOP_ADD)
4113 mpi_add(mpr_num(bigr), &iop, mpr_num(bigr));
4114 else
4115 mpi_sub(mpr_num(bigr), &iop, mpr_num(bigr));
4116 mpi_clear(&iop);
4117 RBR(real) = bigr;
4118 RTYPE(real) = N_BIGRATIO;
4119 rbr_canonicalize(real);
4120 }
4121 }
4122
4123 static void
rop_fi_fr_md_xr(n_real * real,long num,long den,int nop)4124 rop_fi_fr_md_xr(n_real *real, long num, long den, int nop)
4125 {
4126 int fit;
4127 long op = RFI(real);
4128
4129 if (nop == NOP_MUL)
4130 fit = !fi_fi_mul_overflow(op, num);
4131 else
4132 fit = !fi_fi_mul_overflow(op, den);
4133 if (fit) {
4134 if (nop == NOP_MUL) {
4135 RFRN(real) = op * num;
4136 RFRD(real) = den;
4137 }
4138 else {
4139 RFRN(real) = op * den;
4140 RFRD(real) = num;
4141 }
4142 RTYPE(real) = N_FIXRATIO;
4143 rfr_canonicalize(real);
4144 }
4145 else {
4146 mpi iop;
4147 mpr *bigr = XALLOC(mpr);
4148
4149 mpi_init(&iop);
4150 mpi_seti(&iop, op);
4151
4152 mpr_init(bigr);
4153 if (nop == NOP_MUL)
4154 mpr_seti(bigr, num, den);
4155 else
4156 mpr_seti(bigr, den, num);
4157 mpi_mul(mpr_num(bigr), mpr_num(bigr), &iop);
4158 mpi_clear(&iop);
4159 RBR(real) = bigr;
4160 RTYPE(real) = N_BIGRATIO;
4161 rbr_canonicalize(real);
4162 }
4163 }
4164
4165 static INLINE void
radd_fi_fr(n_real * real,long num,long den)4166 radd_fi_fr(n_real *real, long num, long den)
4167 {
4168 rop_fi_fr_as_xr(real, num, den, NOP_ADD);
4169 }
4170
4171 static INLINE void
rsub_fi_fr(n_real * real,long num,long den)4172 rsub_fi_fr(n_real *real, long num, long den)
4173 {
4174 rop_fi_fr_as_xr(real, num, den, NOP_SUB);
4175 }
4176
4177 static INLINE void
rmul_fi_fr(n_real * real,long num,long den)4178 rmul_fi_fr(n_real *real, long num, long den)
4179 {
4180 rop_fi_fr_md_xr(real, num, den, NOP_MUL);
4181 }
4182
4183 static INLINE void
rdiv_fi_fr(n_real * real,long num,long den)4184 rdiv_fi_fr(n_real *real, long num, long den)
4185 {
4186 rop_fi_fr_md_xr(real, num, den, NOP_DIV);
4187 }
4188
4189 static INLINE int
cmp_fi_fr(long fi,long num,long den)4190 cmp_fi_fr(long fi, long num, long den)
4191 {
4192 return (cmp_flonum((double)fi, (double)num / (double)den));
4193 }
4194
4195
4196 /************************************************************************
4197 * FIXNUM BIGRATIO
4198 ************************************************************************/
4199 static void
rop_fi_br_as_xr(n_real * real,mpr * ratio,int nop)4200 rop_fi_br_as_xr(n_real *real, mpr *ratio, int nop)
4201 {
4202 mpi iop;
4203 mpr *bigr = XALLOC(mpr);
4204
4205 mpi_init(&iop);
4206 mpi_seti(&iop, RFI(real));
4207
4208 mpr_init(bigr);
4209 mpr_set(bigr, ratio);
4210
4211 mpi_mul(&iop, &iop, mpr_den(ratio));
4212 if (nop == NOP_ADD)
4213 mpi_add(mpr_num(bigr), &iop, mpr_num(bigr));
4214 else
4215 mpi_sub(mpr_num(bigr), &iop, mpr_num(bigr));
4216
4217 mpi_clear(&iop);
4218 RBR(real) = bigr;
4219 RTYPE(real) = N_BIGRATIO;
4220 rbr_canonicalize(real);
4221 }
4222
4223 static void
rop_fi_br_md_xr(n_real * real,mpr * ratio,int nop)4224 rop_fi_br_md_xr(n_real *real, mpr *ratio, int nop)
4225 {
4226 mpi iop;
4227 mpr *bigr = XALLOC(mpr);
4228
4229 mpi_init(&iop);
4230 mpi_seti(&iop, RFI(real));
4231
4232 mpr_init(bigr);
4233 if (nop == NOP_MUL)
4234 mpr_set(bigr, ratio);
4235 else
4236 mpr_inv(bigr, ratio);
4237
4238 mpi_mul(mpr_num(bigr), &iop, mpr_num(bigr));
4239
4240 mpi_clear(&iop);
4241 RBR(real) = bigr;
4242 RTYPE(real) = N_BIGRATIO;
4243 rbr_canonicalize(real);
4244 }
4245
4246 static INLINE void
radd_fi_br(n_real * real,mpr * ratio)4247 radd_fi_br(n_real *real, mpr *ratio)
4248 {
4249 rop_fi_br_as_xr(real, ratio, NOP_ADD);
4250 }
4251
4252 static INLINE void
rsub_fi_br(n_real * real,mpr * ratio)4253 rsub_fi_br(n_real *real, mpr *ratio)
4254 {
4255 rop_fi_br_as_xr(real, ratio, NOP_SUB);
4256 }
4257
4258 static INLINE void
rmul_fi_br(n_real * real,mpr * ratio)4259 rmul_fi_br(n_real *real, mpr *ratio)
4260 {
4261 rop_fi_br_md_xr(real, ratio, NOP_MUL);
4262 }
4263
4264 static INLINE void
rdiv_fi_br(n_real * real,mpr * ratio)4265 rdiv_fi_br(n_real *real, mpr *ratio)
4266 {
4267 rop_fi_br_md_xr(real, ratio, NOP_DIV);
4268 }
4269
4270 static INLINE int
cmp_fi_br(long op1,mpr * op2)4271 cmp_fi_br(long op1, mpr *op2)
4272 {
4273 return (-mpr_cmpi(op2, op1));
4274 }
4275
4276
4277 /************************************************************************
4278 * BIGNUM FIXNUM
4279 ************************************************************************/
4280 static INLINE void
radd_bi_fi(n_real * real,long fi)4281 radd_bi_fi(n_real *real, long fi)
4282 {
4283 mpi_addi(RBI(real), RBI(real), fi);
4284 rbi_canonicalize(real);
4285 }
4286
4287 static INLINE void
rsub_bi_fi(n_real * real,long fi)4288 rsub_bi_fi(n_real *real, long fi)
4289 {
4290 mpi_subi(RBI(real), RBI(real), fi);
4291 rbi_canonicalize(real);
4292 }
4293
4294 static INLINE void
rmul_bi_fi(n_real * real,long fi)4295 rmul_bi_fi(n_real *real, long fi)
4296 {
4297 mpi_muli(RBI(real), RBI(real), fi);
4298 rbi_canonicalize(real);
4299 }
4300
4301 static void
rdiv_bi_fi(n_real * real,long fi)4302 rdiv_bi_fi(n_real *real, long fi)
4303 {
4304 mpr *bigr;
4305
4306 if (RFI(real) == 0)
4307 fatal_error(DIVIDE_BY_ZERO);
4308
4309 bigr = XALLOC(mpr);
4310 mpr_init(bigr);
4311 mpi_set(mpr_num(bigr), RBI(real));
4312 mpi_seti(mpr_den(bigr), fi);
4313 RCLEAR_BI(real);
4314 RBR(real) = bigr;
4315 RTYPE(real) = N_BIGRATIO;
4316 rbr_canonicalize(real);
4317 }
4318
4319 static INLINE int
cmp_bi_fi(mpi * bignum,long fi)4320 cmp_bi_fi(mpi *bignum, long fi)
4321 {
4322 return (mpi_cmpi(bignum, fi));
4323 }
4324
4325
4326 /************************************************************************
4327 * BIGNUM BIGNUM
4328 ************************************************************************/
4329 static INLINE void
radd_bi_bi(n_real * real,mpi * bignum)4330 radd_bi_bi(n_real *real, mpi *bignum)
4331 {
4332 mpi_add(RBI(real), RBI(real), bignum);
4333 rbi_canonicalize(real);
4334 }
4335
4336 static INLINE void
rsub_bi_bi(n_real * real,mpi * bignum)4337 rsub_bi_bi(n_real *real, mpi *bignum)
4338 {
4339 mpi_sub(RBI(real), RBI(real), bignum);
4340 rbi_canonicalize(real);
4341 }
4342
4343 static INLINE void
rmul_bi_bi(n_real * real,mpi * bignum)4344 rmul_bi_bi(n_real *real, mpi *bignum)
4345 {
4346 mpi_mul(RBI(real), RBI(real), bignum);
4347 rbi_canonicalize(real);
4348 }
4349
4350 static void
rdiv_bi_bi(n_real * real,mpi * bignum)4351 rdiv_bi_bi(n_real *real, mpi *bignum)
4352 {
4353 mpr *bigr;
4354
4355 if (mpi_cmpi(bignum, 0) == 0)
4356 fatal_error(DIVIDE_BY_ZERO);
4357
4358 bigr = XALLOC(mpr);
4359 mpr_init(bigr);
4360 mpi_set(mpr_num(bigr), RBI(real));
4361 mpi_set(mpr_den(bigr), bignum);
4362 RCLEAR_BI(real);
4363 RBR(real) = bigr;
4364 RTYPE(real) = N_BIGRATIO;
4365 rbr_canonicalize(real);
4366 }
4367
4368 static INLINE int
cmp_bi_bi(mpi * op1,mpi * op2)4369 cmp_bi_bi(mpi *op1, mpi *op2)
4370 {
4371 return (mpi_cmp(op1, op2));
4372 }
4373
4374
4375 /************************************************************************
4376 * BIGNUM FIXRATIO
4377 ************************************************************************/
4378 static void
rop_bi_fr_as_xr(n_real * real,long num,long den,int nop)4379 rop_bi_fr_as_xr(n_real *real, long num, long den, int nop)
4380 {
4381 mpi iop;
4382 mpr *bigr = XALLOC(mpr);
4383
4384 mpi_init(&iop);
4385 mpi_set(&iop, RBI(real));
4386 mpi_muli(&iop, &iop, den);
4387
4388 mpr_init(bigr);
4389 mpr_seti(bigr, num, den);
4390
4391 if (nop == NOP_ADD)
4392 mpi_add(mpr_num(bigr), &iop, mpr_num(bigr));
4393 else
4394 mpi_sub(mpr_num(bigr), &iop, mpr_num(bigr));
4395 mpi_clear(&iop);
4396
4397 RCLEAR_BI(real);
4398 RBR(real) = bigr;
4399 RTYPE(real) = N_BIGRATIO;
4400 rbr_canonicalize(real);
4401 }
4402
4403 static INLINE void
rop_bi_fr_md_xr(n_real * real,long num,long den,int nop)4404 rop_bi_fr_md_xr(n_real *real, long num, long den, int nop)
4405 {
4406 mpr *bigr = XALLOC(mpr);
4407
4408 mpr_init(bigr);
4409
4410 mpr_seti(bigr, num, den);
4411
4412 if (nop == NOP_MUL)
4413 mpi_mul(mpr_num(bigr), RBI(real), mpr_num(bigr));
4414 else {
4415 mpi_mul(mpr_den(bigr), RBI(real), mpr_den(bigr));
4416 mpr_inv(bigr, bigr);
4417 }
4418
4419 RCLEAR_BI(real);
4420 RBR(real) = bigr;
4421 RTYPE(real) = N_BIGRATIO;
4422 rbr_canonicalize(real);
4423 }
4424
4425 static INLINE void
radd_bi_fr(n_real * real,long num,long den)4426 radd_bi_fr(n_real *real, long num, long den)
4427 {
4428 rop_bi_fr_as_xr(real, num, den, NOP_ADD);
4429 }
4430
4431 static INLINE void
rsub_bi_fr(n_real * real,long num,long den)4432 rsub_bi_fr(n_real *real, long num, long den)
4433 {
4434 rop_bi_fr_as_xr(real, num, den, NOP_SUB);
4435 }
4436
4437 static INLINE void
rmul_bi_fr(n_real * real,long num,long den)4438 rmul_bi_fr(n_real *real, long num, long den)
4439 {
4440 rop_bi_fr_md_xr(real, num, den, NOP_MUL);
4441 }
4442
4443 static INLINE void
rdiv_bi_fr(n_real * real,long num,long den)4444 rdiv_bi_fr(n_real *real, long num, long den)
4445 {
4446 rop_bi_fr_md_xr(real, num, den, NOP_DIV);
4447 }
4448
4449 static int
cmp_bi_fr(mpi * bignum,long num,long den)4450 cmp_bi_fr(mpi *bignum, long num, long den)
4451 {
4452 int cmp;
4453 mpr cmp1, cmp2;
4454
4455 mpr_init(&cmp1);
4456 mpi_set(mpr_num(&cmp1), bignum);
4457 mpi_seti(mpr_den(&cmp1), 1);
4458
4459 mpr_init(&cmp2);
4460 mpr_seti(&cmp2, num, den);
4461
4462 cmp = mpr_cmp(&cmp1, &cmp2);
4463 mpr_clear(&cmp1);
4464 mpr_clear(&cmp2);
4465
4466 return (cmp);
4467 }
4468
4469
4470 /************************************************************************
4471 * BIGNUM BIGRATIO
4472 ************************************************************************/
4473 static void
rop_bi_br_as_xr(n_real * real,mpr * bigratio,int nop)4474 rop_bi_br_as_xr(n_real *real, mpr *bigratio, int nop)
4475 {
4476 mpi iop;
4477 mpr *bigr = XALLOC(mpr);
4478
4479 mpi_init(&iop);
4480 mpi_set(&iop, RBI(real));
4481 mpr_init(bigr);
4482 mpr_set(bigr, bigratio);
4483
4484 mpi_mul(&iop, &iop, mpr_den(bigratio));
4485
4486 if (nop == NOP_ADD)
4487 mpi_add(mpr_num(bigr), &iop, mpr_num(bigr));
4488 else
4489 mpi_sub(mpr_num(bigr), &iop, mpr_num(bigr));
4490 mpi_clear(&iop);
4491
4492 RCLEAR_BI(real);
4493 RBR(real) = bigr;
4494 RTYPE(real) = N_BIGRATIO;
4495 rbr_canonicalize(real);
4496 }
4497
4498 static void
rop_bi_br_md_xr(n_real * real,mpr * bigratio,int nop)4499 rop_bi_br_md_xr(n_real *real, mpr *bigratio, int nop)
4500 {
4501 mpr *bigr = XALLOC(mpr);
4502
4503 mpr_init(bigr);
4504 if (nop == NOP_MUL)
4505 mpr_set(bigr, bigratio);
4506 else
4507 mpr_inv(bigr, bigratio);
4508
4509 mpi_mul(mpr_num(bigr), RBI(real), mpr_num(bigr));
4510
4511 RCLEAR_BI(real);
4512 RBR(real) = bigr;
4513 RTYPE(real) = N_BIGRATIO;
4514 rbr_canonicalize(real);
4515 }
4516
4517 static INLINE void
radd_bi_br(n_real * real,mpr * bigratio)4518 radd_bi_br(n_real *real, mpr *bigratio)
4519 {
4520 rop_bi_br_as_xr(real, bigratio, NOP_ADD);
4521 }
4522
4523 static INLINE void
rsub_bi_br(n_real * real,mpr * bigratio)4524 rsub_bi_br(n_real *real, mpr *bigratio)
4525 {
4526 rop_bi_br_as_xr(real, bigratio, NOP_SUB);
4527 }
4528
4529 static INLINE void
rmul_bi_br(n_real * real,mpr * bigratio)4530 rmul_bi_br(n_real *real, mpr *bigratio)
4531 {
4532 rop_bi_br_md_xr(real, bigratio, NOP_MUL);
4533 }
4534
4535 static INLINE void
rdiv_bi_br(n_real * real,mpr * bigratio)4536 rdiv_bi_br(n_real *real, mpr *bigratio)
4537 {
4538 rop_bi_br_md_xr(real, bigratio, NOP_DIV);
4539 }
4540
4541 static int
cmp_bi_br(mpi * bignum,mpr * bigratio)4542 cmp_bi_br(mpi *bignum, mpr *bigratio)
4543 {
4544 int cmp;
4545 mpr cmp1;
4546
4547 mpr_init(&cmp1);
4548 mpi_set(mpr_num(&cmp1), bignum);
4549 mpi_seti(mpr_den(&cmp1), 1);
4550
4551 cmp = mpr_cmp(&cmp1, bigratio);
4552 mpr_clear(&cmp1);
4553
4554 return (cmp);
4555 }
4556
4557
4558 /************************************************************************
4559 * FIXRATIO FIXNUM
4560 ************************************************************************/
4561 static void
rop_fr_fi_as_xr(n_real * real,long op,int nop)4562 rop_fr_fi_as_xr(n_real *real, long op, int nop)
4563 {
4564 int fit;
4565 long value = 0, num = RFRN(real), den = RFRD(real);
4566
4567 fit = !fi_fi_mul_overflow(op, den);
4568
4569 if (fit) {
4570 value = op * den;
4571 if (nop == NOP_ADD)
4572 fit = !fi_fi_add_overflow(value, num);
4573 else
4574 fit = !fi_fi_sub_overflow(value, num);
4575 }
4576 if (fit) {
4577 if (nop == NOP_ADD)
4578 RFRN(real) = num + value;
4579 else
4580 RFRN(real) = num - value;
4581 rfr_canonicalize(real);
4582 }
4583 else {
4584 mpi iop;
4585 mpr *bigr = XALLOC(mpr);
4586
4587 mpr_init(bigr);
4588 mpr_seti(bigr, num, den);
4589 mpi_init(&iop);
4590 mpi_seti(&iop, op);
4591 mpi_muli(&iop, &iop, den);
4592 if (nop == NOP_ADD)
4593 mpi_add(mpr_num(bigr), mpr_num(bigr), &iop);
4594 else
4595 mpi_sub(mpr_num(bigr), mpr_num(bigr), &iop);
4596 mpi_clear(&iop);
4597 RBR(real) = bigr;
4598 RTYPE(real) = N_BIGRATIO;
4599 rbr_canonicalize(real);
4600 }
4601 }
4602
4603 static void
rop_fr_fi_md_xr(n_real * real,long op,int nop)4604 rop_fr_fi_md_xr(n_real *real, long op, int nop)
4605 {
4606 long num = RFRN(real), den = RFRD(real);
4607
4608 if (nop == NOP_MUL) {
4609 if (!fi_fi_mul_overflow(op, num)) {
4610 RFRN(real) = op * num;
4611 rfr_canonicalize(real);
4612 return;
4613 }
4614 }
4615 else if (!fi_fi_mul_overflow(op, den)) {
4616 RFRD(real) = op * den;
4617 rfr_canonicalize(real);
4618 return;
4619 }
4620
4621 {
4622 mpr *bigr = XALLOC(mpr);
4623
4624 mpr_init(bigr);
4625 mpr_seti(bigr, num, den);
4626 if (nop == NOP_MUL)
4627 mpr_muli(bigr, bigr, op);
4628 else
4629 mpr_divi(bigr, bigr, op);
4630 RBR(real) = bigr;
4631 RTYPE(real) = N_BIGRATIO;
4632 rbr_canonicalize(real);
4633 }
4634 }
4635
4636 static INLINE void
radd_fr_fi(n_real * real,long op)4637 radd_fr_fi(n_real *real, long op)
4638 {
4639 rop_fr_fi_as_xr(real, op, NOP_ADD);
4640 }
4641
4642 static INLINE void
rsub_fr_fi(n_real * real,long op)4643 rsub_fr_fi(n_real *real, long op)
4644 {
4645 rop_fr_fi_as_xr(real, op, NOP_SUB);
4646 }
4647
4648 static INLINE void
rmul_fr_fi(n_real * real,long op)4649 rmul_fr_fi(n_real *real, long op)
4650 {
4651 rop_fr_fi_md_xr(real, op, NOP_MUL);
4652 }
4653
4654 static INLINE void
rdiv_fr_fi(n_real * real,long op)4655 rdiv_fr_fi(n_real *real, long op)
4656 {
4657 rop_fr_fi_md_xr(real, op, NOP_DIV);
4658 }
4659
4660 static INLINE int
cmp_fr_fi(long num,long den,long fixnum)4661 cmp_fr_fi(long num, long den, long fixnum)
4662 {
4663 return (cmp_flonum((double)num / (double)den, (double)fixnum));
4664 }
4665
4666
4667 /************************************************************************
4668 * FIXRATIO BIGNUM
4669 ************************************************************************/
4670 static void
rop_fr_bi_as_xr(n_real * real,mpi * bignum,int nop)4671 rop_fr_bi_as_xr(n_real *real, mpi *bignum, int nop)
4672 {
4673 mpi iop;
4674 mpr *bigr = XALLOC(mpr);
4675
4676 mpr_init(bigr);
4677 mpr_seti(bigr, RFRN(real), RFRD(real));
4678
4679 mpi_init(&iop);
4680 mpi_set(&iop, bignum);
4681 mpi_muli(&iop, &iop, RFRD(real));
4682
4683 if (nop == NOP_ADD)
4684 mpi_add(mpr_num(bigr), mpr_num(bigr), &iop);
4685 else
4686 mpi_sub(mpr_num(bigr), mpr_num(bigr), &iop);
4687 mpi_clear(&iop);
4688
4689 RBR(real) = bigr;
4690 RTYPE(real) = N_BIGRATIO;
4691 rbr_canonicalize(real);
4692 }
4693
4694 static void
rop_fr_bi_md_xr(n_real * real,mpi * bignum,int nop)4695 rop_fr_bi_md_xr(n_real *real, mpi *bignum, int nop)
4696 {
4697 mpr *bigr = XALLOC(mpr);
4698
4699 mpr_init(bigr);
4700 mpr_seti(bigr, RFRN(real), RFRD(real));
4701
4702 if (nop == NOP_MUL)
4703 mpi_mul(mpr_num(bigr), mpr_num(bigr), bignum);
4704 else
4705 mpi_mul(mpr_den(bigr), mpr_den(bigr), bignum);
4706
4707 RBR(real) = bigr;
4708 RTYPE(real) = N_BIGRATIO;
4709 rbr_canonicalize(real);
4710 }
4711
4712 static INLINE void
radd_fr_bi(n_real * real,mpi * bignum)4713 radd_fr_bi(n_real *real, mpi *bignum)
4714 {
4715 rop_fr_bi_as_xr(real, bignum, NOP_ADD);
4716 }
4717
4718 static INLINE void
rsub_fr_bi(n_real * real,mpi * bignum)4719 rsub_fr_bi(n_real *real, mpi *bignum)
4720 {
4721 rop_fr_bi_as_xr(real, bignum, NOP_SUB);
4722 }
4723
4724 static INLINE void
rmul_fr_bi(n_real * real,mpi * bignum)4725 rmul_fr_bi(n_real *real, mpi *bignum)
4726 {
4727 rop_fr_bi_md_xr(real, bignum, NOP_MUL);
4728 }
4729
4730 static INLINE void
rdiv_fr_bi(n_real * real,mpi * bignum)4731 rdiv_fr_bi(n_real *real, mpi *bignum)
4732 {
4733 rop_fr_bi_md_xr(real, bignum, NOP_DIV);
4734 }
4735
4736 static int
cmp_fr_bi(long num,long den,mpi * bignum)4737 cmp_fr_bi(long num, long den, mpi *bignum)
4738 {
4739 int cmp;
4740 mpr cmp1, cmp2;
4741
4742 mpr_init(&cmp1);
4743 mpr_seti(&cmp1, num, den);
4744
4745 mpr_init(&cmp2);
4746 mpi_set(mpr_num(&cmp2), bignum);
4747 mpi_seti(mpr_den(&cmp2), 1);
4748
4749 cmp = mpr_cmp(&cmp1, &cmp2);
4750 mpr_clear(&cmp1);
4751 mpr_clear(&cmp2);
4752
4753 return (cmp);
4754 }
4755
4756
4757 /************************************************************************
4758 * FIXRATIO FIXRATIO
4759 ************************************************************************/
4760 static void
rop_fr_fr_as_xr(n_real * real,long num2,long den2,int nop)4761 rop_fr_fr_as_xr(n_real *real, long num2, long den2, int nop)
4762 {
4763 int fit;
4764 long num1 = RFRN(real), den1 = RFRD(real), num = 0, den = 0;
4765
4766 fit = !fi_fi_mul_overflow(num1, den2);
4767 if (fit) {
4768 num = num1 * den2;
4769 fit = !fi_fi_mul_overflow(num2, den1);
4770 if (fit) {
4771 den = num2 * den1;
4772 if (nop == NOP_ADD) {
4773 if ((fit = !fi_fi_add_overflow(num, den)) != 0)
4774 num += den;
4775 }
4776 else if ((fit = !fi_fi_sub_overflow(num, den)) != 0)
4777 num -= den;
4778 if (fit) {
4779 fit = !fi_fi_mul_overflow(den1, den2);
4780 if (fit)
4781 den = den1 * den2;
4782 }
4783 }
4784 }
4785 if (fit) {
4786 RFRN(real) = num;
4787 RFRD(real) = den;
4788 rfr_canonicalize(real);
4789 }
4790 else {
4791 mpi iop;
4792 mpr *bigr = XALLOC(mpr);
4793
4794 mpr_init(bigr);
4795 mpr_seti(bigr, num1, den1);
4796 mpi_muli(mpr_den(bigr), mpr_den(bigr), den2);
4797 mpi_init(&iop);
4798 mpi_seti(&iop, num2);
4799 mpi_muli(&iop, &iop, den1);
4800 mpi_muli(mpr_num(bigr), mpr_num(bigr), den2);
4801 if (nop == NOP_ADD)
4802 mpi_add(mpr_num(bigr), mpr_num(bigr), &iop);
4803 else
4804 mpi_sub(mpr_num(bigr), mpr_num(bigr), &iop);
4805 mpi_clear(&iop);
4806 RBR(real) = bigr;
4807 RTYPE(real) = N_BIGRATIO;
4808 rbr_canonicalize(real);
4809 }
4810 }
4811
4812 static void
rop_fr_fr_md_xr(n_real * real,long num2,long den2,int nop)4813 rop_fr_fr_md_xr(n_real *real, long num2, long den2, int nop)
4814 {
4815 int fit;
4816 long num1 = RFRN(real), den1 = RFRD(real), num = 0, den = 0;
4817
4818 if (nop == NOP_MUL) {
4819 fit = !fi_fi_mul_overflow(num1, num2) && !fi_fi_mul_overflow(den1, den2);
4820 if (fit) {
4821 num = num1 * num2;
4822 den = den1 * den2;
4823 }
4824 }
4825 else {
4826 fit = !fi_fi_mul_overflow(num1, den2) && !fi_fi_mul_overflow(den1, num2);
4827 if (fit) {
4828 num = num1 * den2;
4829 den = den1 * num2;
4830 }
4831 }
4832
4833 if (fit) {
4834 RFRN(real) = num;
4835 RFRD(real) = den;
4836 rfr_canonicalize(real);
4837 }
4838 else {
4839 mpr *bigr = XALLOC(mpr);
4840
4841 mpr_init(bigr);
4842
4843 if (nop == NOP_MUL) {
4844 mpr_seti(bigr, num1, den1);
4845 mpi_muli(mpr_num(bigr), mpr_num(bigr), num2);
4846 mpi_muli(mpr_den(bigr), mpr_den(bigr), den2);
4847 }
4848 else {
4849 mpr_seti(bigr, num1, num2);
4850 mpi_muli(mpr_num(bigr), mpr_num(bigr), den2);
4851 mpi_muli(mpr_den(bigr), mpr_den(bigr), den1);
4852 }
4853
4854 RBR(real) = bigr;
4855 RTYPE(real) = N_BIGRATIO;
4856 rbr_canonicalize(real);
4857 }
4858 }
4859
4860 static INLINE void
radd_fr_fr(n_real * real,long num,long den)4861 radd_fr_fr(n_real *real, long num, long den)
4862 {
4863 rop_fr_fr_as_xr(real, num, den, NOP_ADD);
4864 }
4865
4866 static INLINE void
rsub_fr_fr(n_real * real,long num,long den)4867 rsub_fr_fr(n_real *real, long num, long den)
4868 {
4869 rop_fr_fr_as_xr(real, num, den, NOP_SUB);
4870 }
4871
4872 static INLINE void
rmul_fr_fr(n_real * real,long num,long den)4873 rmul_fr_fr(n_real *real, long num, long den)
4874 {
4875 rop_fr_fr_md_xr(real, num, den, NOP_MUL);
4876 }
4877
4878 static INLINE void
rdiv_fr_fr(n_real * real,long num,long den)4879 rdiv_fr_fr(n_real *real, long num, long den)
4880 {
4881 rop_fr_fr_md_xr(real, num, den, NOP_DIV);
4882 }
4883
4884 static INLINE int
cmp_fr_fr(long num1,long den1,long num2,long den2)4885 cmp_fr_fr(long num1, long den1, long num2, long den2)
4886 {
4887 return (cmp_flonum((double)num1 / (double)den1,
4888 (double)num2 / (double)den2));
4889 }
4890
4891
4892 /************************************************************************
4893 * FIXRATIO BIGRATIO
4894 ************************************************************************/
4895 static void
rop_fr_br_asmd_xr(n_real * real,mpr * bigratio,int nop)4896 rop_fr_br_asmd_xr(n_real *real, mpr *bigratio, int nop)
4897 {
4898 mpr *bigr = XALLOC(mpr);
4899
4900 mpr_init(bigr);
4901 mpr_seti(bigr, RFRN(real), RFRD(real));
4902
4903 switch (nop) {
4904 case NOP_ADD:
4905 mpr_add(bigr, bigr, bigratio);
4906 break;
4907 case NOP_SUB:
4908 mpr_sub(bigr, bigr, bigratio);
4909 break;
4910 case NOP_MUL:
4911 mpr_mul(bigr, bigr, bigratio);
4912 break;
4913 default:
4914 mpr_div(bigr, bigr, bigratio);
4915 break;
4916 }
4917
4918 RBR(real) = bigr;
4919 RTYPE(real) = N_BIGRATIO;
4920 rbr_canonicalize(real);
4921 }
4922
4923 static INLINE void
radd_fr_br(n_real * real,mpr * bigratio)4924 radd_fr_br(n_real *real, mpr *bigratio)
4925 {
4926 rop_fr_br_asmd_xr(real, bigratio, NOP_ADD);
4927 }
4928
4929 static INLINE void
rsub_fr_br(n_real * real,mpr * bigratio)4930 rsub_fr_br(n_real *real, mpr *bigratio)
4931 {
4932 rop_fr_br_asmd_xr(real, bigratio, NOP_SUB);
4933 }
4934
4935 static INLINE void
rmul_fr_br(n_real * real,mpr * bigratio)4936 rmul_fr_br(n_real *real, mpr *bigratio)
4937 {
4938 rop_fr_br_asmd_xr(real, bigratio, NOP_MUL);
4939 }
4940
4941 static INLINE void
rdiv_fr_br(n_real * real,mpr * bigratio)4942 rdiv_fr_br(n_real *real, mpr *bigratio)
4943 {
4944 rop_fr_br_asmd_xr(real, bigratio, NOP_DIV);
4945 }
4946
4947 static int
cmp_fr_br(long num,long den,mpr * bigratio)4948 cmp_fr_br(long num, long den, mpr *bigratio)
4949 {
4950 int cmp;
4951 mpr cmp1;
4952
4953 mpr_init(&cmp1);
4954 mpr_seti(&cmp1, num, den);
4955
4956 cmp = mpr_cmp(&cmp1, bigratio);
4957 mpr_clear(&cmp1);
4958
4959 return (cmp);
4960 }
4961
4962
4963 /************************************************************************
4964 * BIGRATIO FIXNUM
4965 ************************************************************************/
4966 static void
rop_br_fi_asmd_xr(n_real * real,long fixnum,int nop)4967 rop_br_fi_asmd_xr(n_real *real, long fixnum, int nop)
4968 {
4969 mpr *bigratio = RBR(real);
4970
4971 switch (nop) {
4972 case NOP_ADD:
4973 mpr_addi(bigratio, bigratio, fixnum);
4974 break;
4975 case NOP_SUB:
4976 mpr_subi(bigratio, bigratio, fixnum);
4977 break;
4978 case NOP_MUL:
4979 mpr_muli(bigratio, bigratio, fixnum);
4980 break;
4981 default:
4982 if (fixnum == 0)
4983 fatal_error(DIVIDE_BY_ZERO);
4984 mpr_divi(bigratio, bigratio, fixnum);
4985 break;
4986 }
4987 rbr_canonicalize(real);
4988 }
4989
4990 static INLINE void
radd_br_fi(n_real * real,long fixnum)4991 radd_br_fi(n_real *real, long fixnum)
4992 {
4993 rop_br_fi_asmd_xr(real, fixnum, NOP_ADD);
4994 }
4995
4996 static INLINE void
rsub_br_fi(n_real * real,long fixnum)4997 rsub_br_fi(n_real *real, long fixnum)
4998 {
4999 rop_br_fi_asmd_xr(real, fixnum, NOP_SUB);
5000 }
5001
5002 static INLINE void
rmul_br_fi(n_real * real,long fixnum)5003 rmul_br_fi(n_real *real, long fixnum)
5004 {
5005 rop_br_fi_asmd_xr(real, fixnum, NOP_MUL);
5006 }
5007
5008 static INLINE void
rdiv_br_fi(n_real * real,long fixnum)5009 rdiv_br_fi(n_real *real, long fixnum)
5010 {
5011 rop_br_fi_asmd_xr(real, fixnum, NOP_DIV);
5012 }
5013
5014 static int
cmp_br_fi(mpr * bigratio,long fixnum)5015 cmp_br_fi(mpr *bigratio, long fixnum)
5016 {
5017 int cmp;
5018 mpr cmp2;
5019
5020 mpr_init(&cmp2);
5021 mpr_seti(&cmp2, fixnum, 1);
5022 cmp = mpr_cmp(bigratio, &cmp2);
5023 mpr_clear(&cmp2);
5024
5025 return (cmp);
5026 }
5027
5028
5029 /************************************************************************
5030 * BIGRATIO BIGNUM
5031 ************************************************************************/
5032 static void
rop_br_bi_as_xr(n_real * real,mpi * bignum,int nop)5033 rop_br_bi_as_xr(n_real *real, mpi *bignum, int nop)
5034 {
5035 mpi iop;
5036
5037 mpi_init(&iop);
5038 mpi_set(&iop, bignum);
5039
5040 mpi_mul(&iop, &iop, RBRD(real));
5041 if (nop == NOP_ADD)
5042 mpi_add(RBRN(real), RBRN(real), &iop);
5043 else
5044 mpi_sub(RBRN(real), RBRN(real), &iop);
5045 mpi_clear(&iop);
5046 rbr_canonicalize(real);
5047 }
5048
5049 static INLINE void
radd_br_bi(n_real * real,mpi * bignum)5050 radd_br_bi(n_real *real, mpi *bignum)
5051 {
5052 rop_br_bi_as_xr(real, bignum, NOP_ADD);
5053 }
5054
5055 static INLINE void
rsub_br_bi(n_real * real,mpi * bignum)5056 rsub_br_bi(n_real *real, mpi *bignum)
5057 {
5058 rop_br_bi_as_xr(real, bignum, NOP_SUB);
5059 }
5060
5061 static INLINE void
rmul_br_bi(n_real * real,mpi * bignum)5062 rmul_br_bi(n_real *real, mpi *bignum)
5063 {
5064 mpi_mul(RBRN(real), RBRN(real), bignum);
5065 rbr_canonicalize(real);
5066 }
5067
5068 static INLINE void
rdiv_br_bi(n_real * real,mpi * bignum)5069 rdiv_br_bi(n_real *real, mpi *bignum)
5070 {
5071 mpi_mul(RBRD(real), RBRD(real), bignum);
5072 rbr_canonicalize(real);
5073 }
5074
5075 static int
cmp_br_bi(mpr * bigratio,mpi * bignum)5076 cmp_br_bi(mpr *bigratio, mpi *bignum)
5077 {
5078 int cmp;
5079 mpr cmp1;
5080
5081 mpr_init(&cmp1);
5082 mpi_set(mpr_num(&cmp1), bignum);
5083 mpi_seti(mpr_den(&cmp1), 1);
5084
5085 cmp = mpr_cmp(bigratio, &cmp1);
5086 mpr_clear(&cmp1);
5087
5088 return (cmp);
5089 }
5090
5091
5092 /************************************************************************
5093 * BIGRATIO FIXRATIO
5094 ************************************************************************/
5095 static void
rop_br_fr_asmd_xr(n_real * real,long num,long den,int nop)5096 rop_br_fr_asmd_xr(n_real *real, long num, long den, int nop)
5097 {
5098 mpr *bigratio = RBR(real), rop;
5099
5100 mpr_init(&rop);
5101 mpr_seti(&rop, num, den);
5102 switch (nop) {
5103 case NOP_ADD:
5104 mpr_add(bigratio, bigratio, &rop);
5105 break;
5106 case NOP_SUB:
5107 mpr_sub(bigratio, bigratio, &rop);
5108 break;
5109 case NOP_MUL:
5110 mpr_mul(bigratio, bigratio, &rop);
5111 break;
5112 default:
5113 mpr_div(bigratio, bigratio, &rop);
5114 break;
5115 }
5116 mpr_clear(&rop);
5117 rbr_canonicalize(real);
5118 }
5119
5120 static INLINE void
radd_br_fr(n_real * real,long num,long den)5121 radd_br_fr(n_real *real, long num, long den)
5122 {
5123 rop_br_fr_asmd_xr(real, num, den, NOP_ADD);
5124 }
5125
5126 static INLINE void
rsub_br_fr(n_real * real,long num,long den)5127 rsub_br_fr(n_real *real, long num, long den)
5128 {
5129 rop_br_fr_asmd_xr(real, num, den, NOP_SUB);
5130 }
5131
5132 static INLINE void
rmul_br_fr(n_real * real,long num,long den)5133 rmul_br_fr(n_real *real, long num, long den)
5134 {
5135 rop_br_fr_asmd_xr(real, num, den, NOP_MUL);
5136 }
5137
5138 static INLINE void
rdiv_br_fr(n_real * real,long num,long den)5139 rdiv_br_fr(n_real *real, long num, long den)
5140 {
5141 rop_br_fr_asmd_xr(real, num, den, NOP_DIV);
5142 }
5143
5144 static int
cmp_br_fr(mpr * bigratio,long num,long den)5145 cmp_br_fr(mpr *bigratio, long num, long den)
5146 {
5147 int cmp;
5148 mpr cmp2;
5149
5150 mpr_init(&cmp2);
5151 mpr_seti(&cmp2, num, den);
5152 cmp = mpr_cmp(bigratio, &cmp2);
5153 mpr_clear(&cmp2);
5154
5155 return (cmp);
5156 }
5157
5158
5159 /************************************************************************
5160 * BIGRATIO BIGRATIO
5161 ************************************************************************/
5162 static INLINE void
radd_br_br(n_real * real,mpr * bigratio)5163 radd_br_br(n_real *real, mpr *bigratio)
5164 {
5165 mpr_add(RBR(real), RBR(real), bigratio);
5166 rbr_canonicalize(real);
5167 }
5168
5169 static INLINE void
rsub_br_br(n_real * real,mpr * bigratio)5170 rsub_br_br(n_real *real, mpr *bigratio)
5171 {
5172 mpr_sub(RBR(real), RBR(real), bigratio);
5173 rbr_canonicalize(real);
5174 }
5175
5176 static INLINE void
rmul_br_br(n_real * real,mpr * bigratio)5177 rmul_br_br(n_real *real, mpr *bigratio)
5178 {
5179 mpr_mul(RBR(real), RBR(real), bigratio);
5180 rbr_canonicalize(real);
5181 }
5182
5183 static INLINE void
rdiv_br_br(n_real * real,mpr * bigratio)5184 rdiv_br_br(n_real *real, mpr *bigratio)
5185 {
5186 mpr_div(RBR(real), RBR(real), bigratio);
5187 rbr_canonicalize(real);
5188 }
5189
5190 static INLINE int
cmp_br_br(mpr * op1,mpr * op2)5191 cmp_br_br(mpr *op1, mpr *op2)
5192 {
5193 return (mpr_cmp(op1, op2));
5194 }
5195