« back to all changes in this revision
Viewing changes to lib/number.c
- browse files at revision 1
- compare with another revision
- download diff
- download tarball
- view history from revision 1
- Committer: Bazaar Package Importer
- Author(s): Dirk Eddelbuettel
- Date: 2002-04-13 11:33:49 UTC
- Revision ID: james.westby@ubuntu.com-20020413113349-hl2r1t730b91ov68
Tags: upstream-1.06
Import upstream version 1.06
- files added:
- Examples
- Test
- bc
- dc
- doc
- h
- lib
added
removed
lib/number.c
1
/* number.c: Implements arbitrary precision numbers. */
2
/*
3
Copyright (C) 1991, 1992, 1993, 1994, 1997, 2000 Free Software Foundation, Inc.
4
5
This program is free software; you can redistribute it and/or modify
6
it under the terms of the GNU General Public License as published by
7
the Free Software Foundation; either version 2 of the License , or
8
(at your option) any later version.
9
10
This program is distributed in the hope that it will be useful,
11
but WITHOUT ANY WARRANTY; without even the implied warranty of
12
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13
GNU General Public License for more details.
14
15
You should have received a copy of the GNU General Public License
16
along with this program; see the file COPYING. If not, write to:
17
18
The Free Software Foundation, Inc.
19
59 Temple Place, Suite 330
20
Boston, MA 02111-1307 USA.
21
22
23
You may contact the author by:
24
e-mail: philnelson@acm.org
25
us-mail: Philip A. Nelson
26
Computer Science Department, 9062
27
Western Washington University
28
Bellingham, WA 98226-9062
29
30
*************************************************************************/
31
32
#include <stdio.h>
33
#include <config.h>
34
#include <number.h>
35
#include <assert.h>
36
#include <stdlib.h>
37
#include <ctype.h>/* Prototypes needed for external utility routines. */
38
39
#define bc_rt_warn rt_warn
40
#define bc_rt_error rt_error
41
#define bc_out_of_memory out_of_memory
42
43
_PROTOTYPE(void rt_warn, (char *mesg ,...));
44
_PROTOTYPE(void rt_error, (char *mesg ,...));
45
_PROTOTYPE(void out_of_memory, (void));
46
47
/* Storage used for special numbers. */
48
bc_num _zero_;
49
bc_num _one_;
50
bc_num _two_;
51
52
static bc_num _bc_Free_list = NULL;
53
54
/* new_num allocates a number and sets fields to known values. */
55
56
bc_num
57
bc_new_num (length, scale)
58
int length, scale;
59
{
60
bc_num temp;
61
62
if (_bc_Free_list != NULL) {
63
temp = _bc_Free_list;
64
_bc_Free_list = temp->n_next;
65
} else {
66
temp = (bc_num) malloc (sizeof(bc_struct));
67
if (temp == NULL) bc_out_of_memory ();
68
}
69
temp->n_sign = PLUS;
70
temp->n_len = length;
71
temp->n_scale = scale;
72
temp->n_refs = 1;
73
temp->n_ptr = (char *) malloc (length+scale);
74
if (temp->n_ptr == NULL) bc_out_of_memory();
75
temp->n_value = temp->n_ptr;
76
memset (temp->n_ptr, 0, length+scale);
77
return temp;
78
}
79
80
/* "Frees" a bc_num NUM. Actually decreases reference count and only
81
frees the storage if reference count is zero. */
82
83
void
84
bc_free_num (num)
85
bc_num *num;
86
{
87
if (*num == NULL) return;
88
(*num)->n_refs--;
89
if ((*num)->n_refs == 0) {
90
if ((*num)->n_ptr)
91
free ((*num)->n_ptr);
92
(*num)->n_next = _bc_Free_list;
93
_bc_Free_list = *num;
94
}
95
*num = NULL;
96
}
97
98
99
/* Intitialize the number package! */
100
101
void
102
bc_init_numbers ()
103
{
104
_zero_ = bc_new_num (1,0);
105
_one_ = bc_new_num (1,0);
106
_one_->n_value[0] = 1;
107
_two_ = bc_new_num (1,0);
108
_two_->n_value[0] = 2;
109
}
110
111
112
/* Make a copy of a number! Just increments the reference count! */
113
114
bc_num
115
bc_copy_num (num)
116
bc_num num;
117
{
118
num->n_refs++;
119
return num;
120
}
121
122
123
/* Initialize a number NUM by making it a copy of zero. */
124
125
void
126
bc_init_num (num)
127
bc_num *num;
128
{
129
*num = bc_copy_num (_zero_);
130
}
131
132
/* For many things, we may have leading zeros in a number NUM.
133
_bc_rm_leading_zeros just moves the data "value" pointer to the
134
correct place and adjusts the length. */
135
136
static void
137
_bc_rm_leading_zeros (num)
138
bc_num num;
139
{
140
/* We can move n_value to point to the first non zero digit! */
141
while (*num->n_value == 0 && num->n_len > 1) {
142
num->n_value++;
143
num->n_len--;
144
}
145
}
146
147
148
/* Compare two bc numbers. Return value is 0 if equal, -1 if N1 is less
149
than N2 and +1 if N1 is greater than N2. If USE_SIGN is false, just
150
compare the magnitudes. */
151
152
static int
153
_bc_do_compare (n1, n2, use_sign, ignore_last)
154
bc_num n1, n2;
155
int use_sign;
156
int ignore_last;
157
{
158
char *n1ptr, *n2ptr;
159
int count;
160
161
/* First, compare signs. */
162
if (use_sign && n1->n_sign != n2->n_sign)
163
{
164
if (n1->n_sign == PLUS)
165
return (1); /* Positive N1 > Negative N2 */
166
else
167
return (-1); /* Negative N1 < Positive N1 */
168
}
169
170
/* Now compare the magnitude. */
171
if (n1->n_len != n2->n_len)
172
{
173
if (n1->n_len > n2->n_len)
174
{
175
/* Magnitude of n1 > n2. */
176
if (!use_sign || n1->n_sign == PLUS)
177
return (1);
178
else
179
return (-1);
180
}
181
else
182
{
183
/* Magnitude of n1 < n2. */
184
if (!use_sign || n1->n_sign == PLUS)
185
return (-1);
186
else
187
return (1);
188
}
189
}
190
191
/* If we get here, they have the same number of integer digits.
192
check the integer part and the equal length part of the fraction. */
193
count = n1->n_len + MIN (n1->n_scale, n2->n_scale);
194
n1ptr = n1->n_value;
195
n2ptr = n2->n_value;
196
197
while ((count > 0) && (*n1ptr == *n2ptr))
198
{
199
n1ptr++;
200
n2ptr++;
201
count--;
202
}
203
if (ignore_last && count == 1 && n1->n_scale == n2->n_scale)
204
return (0);
205
if (count != 0)
206
{
207
if (*n1ptr > *n2ptr)
208
{
209
/* Magnitude of n1 > n2. */
210
if (!use_sign || n1->n_sign == PLUS)
211
return (1);
212
else
213
return (-1);
214
}
215
else
216
{
217
/* Magnitude of n1 < n2. */
218
if (!use_sign || n1->n_sign == PLUS)
219
return (-1);
220
else
221
return (1);
222
}
223
}
224
225
/* They are equal up to the last part of the equal part of the fraction. */
226
if (n1->n_scale != n2->n_scale)
227
{
228
if (n1->n_scale > n2->n_scale)
229
{
230
for (count = n1->n_scale-n2->n_scale; count>0; count--)
231
if (*n1ptr++ != 0)
232
{
233
/* Magnitude of n1 > n2. */
234
if (!use_sign || n1->n_sign == PLUS)
235
return (1);
236
else
237
return (-1);
238
}
239
}
240
else
241
{
242
for (count = n2->n_scale-n1->n_scale; count>0; count--)
243
if (*n2ptr++ != 0)
244
{
245
/* Magnitude of n1 < n2. */
246
if (!use_sign || n1->n_sign == PLUS)
247
return (-1);
248
else
249
return (1);
250
}
251
}
252
}
253
254
/* They must be equal! */
255
return (0);
256
}
257
258
259
/* This is the "user callable" routine to compare numbers N1 and N2. */
260
261
int
262
bc_compare (n1, n2)
263
bc_num n1, n2;
264
{
265
return _bc_do_compare (n1, n2, TRUE, FALSE);
266
}
267
268
/* In some places we need to check if the number is negative. */
269
270
char
271
bc_is_neg (num)
272
bc_num num;
273
{
274
return num->n_sign == MINUS;
275
}
276
277
/* In some places we need to check if the number NUM is zero. */
278
279
char
280
bc_is_zero (num)
281
bc_num num;
282
{
283
int count;
284
char *nptr;
285
286
/* Quick check. */
287
if (num == _zero_) return TRUE;
288
289
/* Initialize */
290
count = num->n_len + num->n_scale;
291
nptr = num->n_value;
292
293
/* The check */
294
while ((count > 0) && (*nptr++ == 0)) count--;
295
296
if (count != 0)
297
return FALSE;
298
else
299
return TRUE;
300
}
301
302
/* In some places we need to check if the number NUM is almost zero.
303
Specifically, all but the last digit is 0 and the last digit is 1.
304
Last digit is defined by scale. */
305
306
char
307
bc_is_near_zero (num, scale)
308
bc_num num;
309
int scale;
310
{
311
int count;
312
char *nptr;
313
314
/* Error checking */
315
if (scale > num->n_scale)
316
scale = num->n_scale;
317
318
/* Initialize */
319
count = num->n_len + scale;
320
nptr = num->n_value;
321
322
/* The check */
323
while ((count > 0) && (*nptr++ == 0)) count--;
324
325
if (count != 0 && (count != 1 || *--nptr != 1))
326
return FALSE;
327
else
328
return TRUE;
329
}
330
331
332
/* Perform addition: N1 is added to N2 and the value is
333
returned. The signs of N1 and N2 are ignored.
334
SCALE_MIN is to set the minimum scale of the result. */
335
336
static bc_num
337
_bc_do_add (n1, n2, scale_min)
338
bc_num n1, n2;
339
int scale_min;
340
{
341
bc_num sum;
342
int sum_scale, sum_digits;
343
char *n1ptr, *n2ptr, *sumptr;
344
int carry, n1bytes, n2bytes;
345
int count;
346
347
/* Prepare sum. */
348
sum_scale = MAX (n1->n_scale, n2->n_scale);
349
sum_digits = MAX (n1->n_len, n2->n_len) + 1;
350
sum = bc_new_num (sum_digits, MAX(sum_scale, scale_min));
351
352
/* Zero extra digits made by scale_min. */
353
if (scale_min > sum_scale)
354
{
355
sumptr = (char *) (sum->n_value + sum_scale + sum_digits);
356
for (count = scale_min - sum_scale; count > 0; count--)
357
*sumptr++ = 0;
358
}
359
360
/* Start with the fraction part. Initialize the pointers. */
361
n1bytes = n1->n_scale;
362
n2bytes = n2->n_scale;
363
n1ptr = (char *) (n1->n_value + n1->n_len + n1bytes - 1);
364
n2ptr = (char *) (n2->n_value + n2->n_len + n2bytes - 1);
365
sumptr = (char *) (sum->n_value + sum_scale + sum_digits - 1);
366
367
/* Add the fraction part. First copy the longer fraction.*/
368
if (n1bytes != n2bytes)
369
{
370
if (n1bytes > n2bytes)
371
while (n1bytes>n2bytes)
372
{ *sumptr-- = *n1ptr--; n1bytes--;}
373
else
374
while (n2bytes>n1bytes)
375
{ *sumptr-- = *n2ptr--; n2bytes--;}
376
}
377
378
/* Now add the remaining fraction part and equal size integer parts. */
379
n1bytes += n1->n_len;
380
n2bytes += n2->n_len;
381
carry = 0;
382
while ((n1bytes > 0) && (n2bytes > 0))
383
{
384
*sumptr = *n1ptr-- + *n2ptr-- + carry;
385
if (*sumptr > (BASE-1))
386
{
387
carry = 1;
388
*sumptr -= BASE;
389
}
390
else
391
carry = 0;
392
sumptr--;
393
n1bytes--;
394
n2bytes--;
395
}
396
397
/* Now add carry the longer integer part. */
398
if (n1bytes == 0)
399
{ n1bytes = n2bytes; n1ptr = n2ptr; }
400
while (n1bytes-- > 0)
401
{
402
*sumptr = *n1ptr-- + carry;
403
if (*sumptr > (BASE-1))
404
{
405
carry = 1;
406
*sumptr -= BASE;
407
}
408
else
409
carry = 0;
410
sumptr--;
411
}
412
413
/* Set final carry. */
414
if (carry == 1)
415
*sumptr += 1;
416
417
/* Adjust sum and return. */
418
_bc_rm_leading_zeros (sum);
419
return sum;
420
}
421
422
423
/* Perform subtraction: N2 is subtracted from N1 and the value is
424
returned. The signs of N1 and N2 are ignored. Also, N1 is
425
assumed to be larger than N2. SCALE_MIN is the minimum scale
426
of the result. */
427
428
static bc_num
429
_bc_do_sub (n1, n2, scale_min)
430
bc_num n1, n2;
431
int scale_min;
432
{
433
bc_num diff;
434
int diff_scale, diff_len;
435
int min_scale, min_len;
436
char *n1ptr, *n2ptr, *diffptr;
437
int borrow, count, val;
438
439
/* Allocate temporary storage. */
440
diff_len = MAX (n1->n_len, n2->n_len);
441
diff_scale = MAX (n1->n_scale, n2->n_scale);
442
min_len = MIN (n1->n_len, n2->n_len);
443
min_scale = MIN (n1->n_scale, n2->n_scale);
444
diff = bc_new_num (diff_len, MAX(diff_scale, scale_min));
445
446
/* Zero extra digits made by scale_min. */
447
if (scale_min > diff_scale)
448
{
449
diffptr = (char *) (diff->n_value + diff_len + diff_scale);
450
for (count = scale_min - diff_scale; count > 0; count--)
451
*diffptr++ = 0;
452
}
453
454
/* Initialize the subtract. */
455
n1ptr = (char *) (n1->n_value + n1->n_len + n1->n_scale -1);
456
n2ptr = (char *) (n2->n_value + n2->n_len + n2->n_scale -1);
457
diffptr = (char *) (diff->n_value + diff_len + diff_scale -1);
458
459
/* Subtract the numbers. */
460
borrow = 0;
461
462
/* Take care of the longer scaled number. */
463
if (n1->n_scale != min_scale)
464
{
465
/* n1 has the longer scale */
466
for (count = n1->n_scale - min_scale; count > 0; count--)
467
*diffptr-- = *n1ptr--;
468
}
469
else
470
{
471
/* n2 has the longer scale */
472
for (count = n2->n_scale - min_scale; count > 0; count--)
473
{
474
val = - *n2ptr-- - borrow;
475
if (val < 0)
476
{
477
val += BASE;
478
borrow = 1;
479
}
480
else
481
borrow = 0;
482
*diffptr-- = val;
483
}
484
}
485
486
/* Now do the equal length scale and integer parts. */
487
488
for (count = 0; count < min_len + min_scale; count++)
489
{
490
val = *n1ptr-- - *n2ptr-- - borrow;
491
if (val < 0)
492
{
493
val += BASE;
494
borrow = 1;
495
}
496
else
497
borrow = 0;
498
*diffptr-- = val;
499
}
500
501
/* If n1 has more digits then n2, we now do that subtract. */
502
if (diff_len != min_len)
503
{
504
for (count = diff_len - min_len; count > 0; count--)
505
{
506
val = *n1ptr-- - borrow;
507
if (val < 0)
508
{
509
val += BASE;
510
borrow = 1;
511
}
512
else
513
borrow = 0;
514
*diffptr-- = val;
515
}
516
}
517
518
/* Clean up and return. */
519
_bc_rm_leading_zeros (diff);
520
return diff;
521
}
522
523
524
/* Here is the full subtract routine that takes care of negative numbers.
525
N2 is subtracted from N1 and the result placed in RESULT. SCALE_MIN
526
is the minimum scale for the result. */
527
528
void
529
bc_sub (n1, n2, result, scale_min)
530
bc_num n1, n2, *result;
531
int scale_min;
532
{
533
bc_num diff = NULL;
534
int cmp_res;
535
int res_scale;
536
537
if (n1->n_sign != n2->n_sign)
538
{
539
diff = _bc_do_add (n1, n2, scale_min);
540
diff->n_sign = n1->n_sign;
541
}
542
else
543
{
544
/* subtraction must be done. */
545
/* Compare magnitudes. */
546
cmp_res = _bc_do_compare (n1, n2, FALSE, FALSE);
547
switch (cmp_res)
548
{
549
case -1:
550
/* n1 is less than n2, subtract n1 from n2. */
551
diff = _bc_do_sub (n2, n1, scale_min);
552
diff->n_sign = (n2->n_sign == PLUS ? MINUS : PLUS);
553
break;
554
case 0:
555
/* They are equal! return zero! */
556
res_scale = MAX (scale_min, MAX(n1->n_scale, n2->n_scale));
557
diff = bc_new_num (1, res_scale);
558
memset (diff->n_value, 0, res_scale+1);
559
break;
560
case 1:
561
/* n2 is less than n1, subtract n2 from n1. */
562
diff = _bc_do_sub (n1, n2, scale_min);
563
diff->n_sign = n1->n_sign;
564
break;
565
}
566
}
567
568
/* Clean up and return. */
569
bc_free_num (result);
570
*result = diff;
571
}
572
573
574
/* Here is the full add routine that takes care of negative numbers.
575
N1 is added to N2 and the result placed into RESULT. SCALE_MIN
576
is the minimum scale for the result. */
577
578
void
579
bc_add (n1, n2, result, scale_min)
580
bc_num n1, n2, *result;
581
int scale_min;
582
{
583
bc_num sum = NULL;
584
int cmp_res;
585
int res_scale;
586
587
if (n1->n_sign == n2->n_sign)
588
{
589
sum = _bc_do_add (n1, n2, scale_min);
590
sum->n_sign = n1->n_sign;
591
}
592
else
593
{
594
/* subtraction must be done. */
595
cmp_res = _bc_do_compare (n1, n2, FALSE, FALSE); /* Compare magnitudes. */
596
switch (cmp_res)
597
{
598
case -1:
599
/* n1 is less than n2, subtract n1 from n2. */
600
sum = _bc_do_sub (n2, n1, scale_min);
601
sum->n_sign = n2->n_sign;
602
break;
603
case 0:
604
/* They are equal! return zero with the correct scale! */
605
res_scale = MAX (scale_min, MAX(n1->n_scale, n2->n_scale));
606
sum = bc_new_num (1, res_scale);
607
memset (sum->n_value, 0, res_scale+1);
608
break;
609
case 1:
610
/* n2 is less than n1, subtract n2 from n1. */
611
sum = _bc_do_sub (n1, n2, scale_min);
612
sum->n_sign = n1->n_sign;
613
}
614
}
615
616
/* Clean up and return. */
617
bc_free_num (result);
618
*result = sum;
619
}
620
621
/* Recursive vs non-recursive multiply crossover ranges. */
622
#if defined(MULDIGITS)
623
#include "muldigits.h"
624
#else
625
#define MUL_BASE_DIGITS 80
626
#endif
627
628
int mul_base_digits = MUL_BASE_DIGITS;
629
#define MUL_SMALL_DIGITS mul_base_digits/4
630
631
/* Multiply utility routines */
632
633
static bc_num
634
new_sub_num (length, scale, value)
635
int length, scale;
636
char *value;
637
{
638
bc_num temp;
639
640
if (_bc_Free_list != NULL) {
641
temp = _bc_Free_list;
642
_bc_Free_list = temp->n_next;
643
} else {
644
temp = (bc_num) malloc (sizeof(bc_struct));
645
if (temp == NULL) bc_out_of_memory ();
646
}
647
temp->n_sign = PLUS;
648
temp->n_len = length;
649
temp->n_scale = scale;
650
temp->n_refs = 1;
651
temp->n_ptr = NULL;
652
temp->n_value = value;
653
return temp;
654
}
655
656
static void
657
_bc_simp_mul (bc_num n1, int n1len, bc_num n2, int n2len, bc_num *prod,
658
int full_scale)
659
{
660
char *n1ptr, *n2ptr, *pvptr;
661
char *n1end, *n2end; /* To the end of n1 and n2. */
662
int indx, sum, prodlen;
663
664
prodlen = n1len+n2len+1;
665
666
*prod = bc_new_num (prodlen, 0);
667
668
n1end = (char *) (n1->n_value + n1len - 1);
669
n2end = (char *) (n2->n_value + n2len - 1);
670
pvptr = (char *) ((*prod)->n_value + prodlen - 1);
671
sum = 0;
672
673
/* Here is the loop... */
674
for (indx = 0; indx < prodlen-1; indx++)
675
{
676
n1ptr = (char *) (n1end - MAX(0, indx-n2len+1));
677
n2ptr = (char *) (n2end - MIN(indx, n2len-1));
678
while ((n1ptr >= n1->n_value) && (n2ptr <= n2end))
679
sum += *n1ptr-- * *n2ptr++;
680
*pvptr-- = sum % BASE;
681
sum = sum / BASE;
682
}
683
*pvptr = sum;
684
}
685
686
687
/* A special adder/subtractor for the recursive divide and conquer
688
multiply algorithm. Note: if sub is called, accum must
689
be larger that what is being subtracted. Also, accum and val
690
must have n_scale = 0. (e.g. they must look like integers. *) */
691
static void
692
_bc_shift_addsub (bc_num accum, bc_num val, int shift, int sub)
693
{
694
signed char *accp, *valp;
695
int count, carry;
696
697
count = val->n_len;
698
if (val->n_value[0] == 0)
699
count--;
700
assert (accum->n_len+accum->n_scale >= shift+count);
701
702
/* Set up pointers and others */
703
accp = (signed char *)(accum->n_value +
704
accum->n_len + accum->n_scale - shift - 1);
705
valp = (signed char *)(val->n_value + val->n_len - 1);
706
carry = 0;
707
708
if (sub) {
709
/* Subtraction, carry is really borrow. */
710
while (count--) {
711
*accp -= *valp-- + carry;
712
if (*accp < 0) {
713
carry = 1;
714
*accp-- += BASE;
715
} else {
716
carry = 0;
717
accp--;
718
}
719
}
720
while (carry) {
721
*accp -= carry;
722
if (*accp < 0)
723
*accp-- += BASE;
724
else
725
carry = 0;
726
}
727
} else {
728
/* Addition */
729
while (count--) {
730
*accp += *valp-- + carry;
731
if (*accp > (BASE-1)) {
732
carry = 1;
733
*accp-- -= BASE;
734
} else {
735
carry = 0;
736
accp--;
737
}
738
}
739
while (carry) {
740
*accp += carry;
741
if (*accp > (BASE-1))
742
*accp-- -= BASE;
743
else
744
carry = 0;
745
}
746
}
747
}
748
749
/* Recursive divide and conquer multiply algorithm.
750
Based on
751
Let u = u0 + u1*(b^n)
752
Let v = v0 + v1*(b^n)
753
Then uv = (B^2n+B^n)*u1*v1 + B^n*(u1-u0)*(v0-v1) + (B^n+1)*u0*v0
754
755
B is the base of storage, number of digits in u1,u0 close to equal.
756
*/
757
static void
758
_bc_rec_mul (bc_num u, int ulen, bc_num v, int vlen, bc_num *prod,
759
int full_scale)
760
{
761
bc_num u0, u1, v0, v1;
762
int u0len, v0len;
763
bc_num m1, m2, m3, d1, d2;
764
int n, prodlen, m1zero;
765
int d1len, d2len;
766
767
/* Base case? */
768
if ((ulen+vlen) < mul_base_digits
769
|| ulen < MUL_SMALL_DIGITS
770
|| vlen < MUL_SMALL_DIGITS ) {
771
_bc_simp_mul (u, ulen, v, vlen, prod, full_scale);
772
return;
773
}
774
775
/* Calculate n -- the u and v split point in digits. */
776
n = (MAX(ulen, vlen)+1) / 2;
777
778
/* Split u and v. */
779
if (ulen < n) {
780
u1 = bc_copy_num (_zero_);
781
u0 = new_sub_num (ulen,0, u->n_value);
782
} else {
783
u1 = new_sub_num (ulen-n, 0, u->n_value);
784
u0 = new_sub_num (n, 0, u->n_value+ulen-n);
785
}
786
if (vlen < n) {
787
v1 = bc_copy_num (_zero_);
788
v0 = new_sub_num (vlen,0, v->n_value);
789
} else {
790
v1 = new_sub_num (vlen-n, 0, v->n_value);
791
v0 = new_sub_num (n, 0, v->n_value+vlen-n);
792
}
793
_bc_rm_leading_zeros (u1);
794
_bc_rm_leading_zeros (u0);
795
u0len = u0->n_len;
796
_bc_rm_leading_zeros (v1);
797
_bc_rm_leading_zeros (v0);
798
v0len = v0->n_len;
799
800
m1zero = bc_is_zero(u1) || bc_is_zero(v1);
801
802
/* Calculate sub results ... */
803
804
bc_init_num(&d1);
805
bc_init_num(&d2);
806
bc_sub (u1, u0, &d1, 0);
807
d1len = d1->n_len;
808
bc_sub (v0, v1, &d2, 0);
809
d2len = d2->n_len;
810
811
812
/* Do recursive multiplies and shifted adds. */
813
if (m1zero)
814
m1 = bc_copy_num (_zero_);
815
else
816
_bc_rec_mul (u1, u1->n_len, v1, v1->n_len, &m1, 0);
817
818
if (bc_is_zero(d1) || bc_is_zero(d2))
819
m2 = bc_copy_num (_zero_);
820
else
821
_bc_rec_mul (d1, d1len, d2, d2len, &m2, 0);
822
823
if (bc_is_zero(u0) || bc_is_zero(v0))
824
m3 = bc_copy_num (_zero_);
825
else
826
_bc_rec_mul (u0, u0->n_len, v0, v0->n_len, &m3, 0);
827
828
/* Initialize product */
829
prodlen = ulen+vlen+1;
830
*prod = bc_new_num(prodlen, 0);
831
832
if (!m1zero) {
833
_bc_shift_addsub (*prod, m1, 2*n, 0);
834
_bc_shift_addsub (*prod, m1, n, 0);
835
}
836
_bc_shift_addsub (*prod, m3, n, 0);
837
_bc_shift_addsub (*prod, m3, 0, 0);
838
_bc_shift_addsub (*prod, m2, n, d1->n_sign != d2->n_sign);
839
840
/* Now clean up! */
841
bc_free_num (&u1);
842
bc_free_num (&u0);
843
bc_free_num (&v1);
844
bc_free_num (&m1);
845
bc_free_num (&v0);
846
bc_free_num (&m2);
847
bc_free_num (&m3);
848
bc_free_num (&d1);
849
bc_free_num (&d2);
850
}
851
852
/* The multiply routine. N2 times N1 is put int PROD with the scale of
853
the result being MIN(N2 scale+N1 scale, MAX (SCALE, N2 scale, N1 scale)).
854
*/
855
856
void
857
bc_multiply (n1, n2, prod, scale)
858
bc_num n1, n2, *prod;
859
int scale;
860
{
861
bc_num pval;
862
int len1, len2;
863
int full_scale, prod_scale;
864
865
/* Initialize things. */
866
len1 = n1->n_len + n1->n_scale;
867
len2 = n2->n_len + n2->n_scale;
868
full_scale = n1->n_scale + n2->n_scale;
869
prod_scale = MIN(full_scale,MAX(scale,MAX(n1->n_scale,n2->n_scale)));
870
871
/* Do the multiply */
872
_bc_rec_mul (n1, len1, n2, len2, &pval, full_scale);
873
874
/* Assign to prod and clean up the number. */
875
pval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS );
876
pval->n_value = pval->n_ptr;
877
pval->n_len = len2 + len1 + 1 - full_scale;
878
pval->n_scale = prod_scale;
879
_bc_rm_leading_zeros (pval);
880
if (bc_is_zero (pval))
881
pval->n_sign = PLUS;
882
bc_free_num (prod);
883
*prod = pval;
884
}
885
886
/* Some utility routines for the divide: First a one digit multiply.
887
NUM (with SIZE digits) is multiplied by DIGIT and the result is
888
placed into RESULT. It is written so that NUM and RESULT can be
889
the same pointers. */
890
891
static void
892
_one_mult (num, size, digit, result)
893
unsigned char *num;
894
int size, digit;
895
unsigned char *result;
896
{
897
int carry, value;
898
unsigned char *nptr, *rptr;
899
900
if (digit == 0)
901
memset (result, 0, size);
902
else
903
{
904
if (digit == 1)
905
memcpy (result, num, size);
906
else
907
{
908
/* Initialize */
909
nptr = (unsigned char *) (num+size-1);
910
rptr = (unsigned char *) (result+size-1);
911
carry = 0;
912
913
while (size-- > 0)
914
{
915
value = *nptr-- * digit + carry;
916
*rptr-- = value % BASE;
917
carry = value / BASE;
918
}
919
920
if (carry != 0) *rptr = carry;
921
}
922
}
923
}
924
925
926
/* The full division routine. This computes N1 / N2. It returns
927
0 if the division is ok and the result is in QUOT. The number of
928
digits after the decimal point is SCALE. It returns -1 if division
929
by zero is tried. The algorithm is found in Knuth Vol 2. p237. */
930
931
int
932
bc_divide (n1, n2, quot, scale)
933
bc_num n1, n2, *quot;
934
int scale;
935
{
936
bc_num qval;
937
unsigned char *num1, *num2;
938
unsigned char *ptr1, *ptr2, *n2ptr, *qptr;
939
int scale1, val;
940
unsigned int len1, len2, scale2, qdigits, extra, count;
941
unsigned int qdig, qguess, borrow, carry;
942
unsigned char *mval;
943
char zero;
944
unsigned int norm;
945
946
/* Test for divide by zero. */
947
if (bc_is_zero (n2)) return -1;
948
949
/* Test for divide by 1. If it is we must truncate. */
950
if (n2->n_scale == 0)
951
{
952
if (n2->n_len == 1 && *n2->n_value == 1)
953
{
954
qval = bc_new_num (n1->n_len, scale);
955
qval->n_sign = (n1->n_sign == n2->n_sign ? PLUS : MINUS);
956
memset (&qval->n_value[n1->n_len],0,scale);
957
memcpy (qval->n_value, n1->n_value,
958
n1->n_len + MIN(n1->n_scale,scale));
959
bc_free_num (quot);
960
*quot = qval;
961
}
962
}
963
964
/* Set up the divide. Move the decimal point on n1 by n2's scale.
965
Remember, zeros on the end of num2 are wasted effort for dividing. */
966
scale2 = n2->n_scale;
967
n2ptr = (unsigned char *) n2->n_value+n2->n_len+scale2-1;
968
while ((scale2 > 0) && (*n2ptr-- == 0)) scale2--;
969
970
len1 = n1->n_len + scale2;
971
scale1 = n1->n_scale - scale2;
972
if (scale1 < scale)
973
extra = scale - scale1;
974
else
975
extra = 0;
976
num1 = (unsigned char *) malloc (n1->n_len+n1->n_scale+extra+2);
977
if (num1 == NULL) bc_out_of_memory();
978
memset (num1, 0, n1->n_len+n1->n_scale+extra+2);
979
memcpy (num1+1, n1->n_value, n1->n_len+n1->n_scale);
980
981
len2 = n2->n_len + scale2;
982
num2 = (unsigned char *) malloc (len2+1);
983
if (num2 == NULL) bc_out_of_memory();
984
memcpy (num2, n2->n_value, len2);
985
*(num2+len2) = 0;
986
n2ptr = num2;
987
while (*n2ptr == 0)
988
{
989
n2ptr++;
990
len2--;
991
}
992
993
/* Calculate the number of quotient digits. */
994
if (len2 > len1+scale)
995
{
996
qdigits = scale+1;
997
zero = TRUE;
998
}
999
else
1000
{
1001
zero = FALSE;
1002
if (len2>len1)
1003
qdigits = scale+1; /* One for the zero integer part. */
1004
else
1005
qdigits = len1-len2+scale+1;
1006
}
1007
1008
/* Allocate and zero the storage for the quotient. */
1009
qval = bc_new_num (qdigits-scale,scale);
1010
memset (qval->n_value, 0, qdigits);
1011
1012
/* Allocate storage for the temporary storage mval. */
1013
mval = (unsigned char *) malloc (len2+1);
1014
if (mval == NULL) bc_out_of_memory ();
1015
1016
/* Now for the full divide algorithm. */
1017
if (!zero)
1018
{
1019
/* Normalize */
1020
norm = 10 / ((int)*n2ptr + 1);
1021
if (norm != 1)
1022
{
1023
_one_mult (num1, len1+scale1+extra+1, norm, num1);
1024
_one_mult (n2ptr, len2, norm, n2ptr);
1025
}
1026
1027
/* Initialize divide loop. */
1028
qdig = 0;
1029
if (len2 > len1)
1030
qptr = (unsigned char *) qval->n_value+len2-len1;
1031
else
1032
qptr = (unsigned char *) qval->n_value;
1033
1034
/* Loop */
1035
while (qdig <= len1+scale-len2)
1036
{
1037
/* Calculate the quotient digit guess. */
1038
if (*n2ptr == num1[qdig])
1039
qguess = 9;
1040
else
1041
qguess = (num1[qdig]*10 + num1[qdig+1]) / *n2ptr;
1042
1043
/* Test qguess. */
1044
if (n2ptr[1]*qguess >
1045
(num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10
1046
+ num1[qdig+2])
1047
{
1048
qguess--;
1049
/* And again. */
1050
if (n2ptr[1]*qguess >
1051
(num1[qdig]*10 + num1[qdig+1] - *n2ptr*qguess)*10
1052
+ num1[qdig+2])
1053
qguess--;
1054
}
1055
1056
/* Multiply and subtract. */
1057
borrow = 0;
1058
if (qguess != 0)
1059
{
1060
*mval = 0;
1061
_one_mult (n2ptr, len2, qguess, mval+1);
1062
ptr1 = (unsigned char *) num1+qdig+len2;
1063
ptr2 = (unsigned char *) mval+len2;
1064
for (count = 0; count < len2+1; count++)
1065
{
1066
val = (int) *ptr1 - (int) *ptr2-- - borrow;
1067
if (val < 0)
1068
{
1069
val += 10;
1070
borrow = 1;
1071
}
1072
else
1073
borrow = 0;
1074
*ptr1-- = val;
1075
}
1076
}
1077
1078
/* Test for negative result. */
1079
if (borrow == 1)
1080
{
1081
qguess--;
1082
ptr1 = (unsigned char *) num1+qdig+len2;
1083
ptr2 = (unsigned char *) n2ptr+len2-1;
1084
carry = 0;
1085
for (count = 0; count < len2; count++)
1086
{
1087
val = (int) *ptr1 + (int) *ptr2-- + carry;
1088
if (val > 9)
1089
{
1090
val -= 10;
1091
carry = 1;
1092
}
1093
else
1094
carry = 0;
1095
*ptr1-- = val;
1096
}
1097
if (carry == 1) *ptr1 = (*ptr1 + 1) % 10;
1098
}
1099
1100
/* We now know the quotient digit. */
1101
*qptr++ = qguess;
1102
qdig++;
1103
}
1104
}
1105
1106
/* Clean up and return the number. */
1107
qval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS );
1108
if (bc_is_zero (qval)) qval->n_sign = PLUS;
1109
_bc_rm_leading_zeros (qval);
1110
bc_free_num (quot);
1111
*quot = qval;
1112
1113
/* Clean up temporary storage. */
1114
free (mval);
1115
free (num1);
1116
free (num2);
1117
1118
return 0; /* Everything is OK. */
1119
}
1120
1121
1122
/* Division *and* modulo for numbers. This computes both NUM1 / NUM2 and
1123
NUM1 % NUM2 and puts the results in QUOT and REM, except that if QUOT
1124
is NULL then that store will be omitted.
1125
*/
1126
1127
int
1128
bc_divmod (num1, num2, quot, rem, scale)
1129
bc_num num1, num2, *quot, *rem;
1130
int scale;
1131
{
1132
bc_num quotient = NULL;
1133
bc_num temp;
1134
int rscale;
1135
1136
/* Check for correct numbers. */
1137
if (bc_is_zero (num2)) return -1;
1138
1139
/* Calculate final scale. */
1140
rscale = MAX (num1->n_scale, num2->n_scale+scale);
1141
bc_init_num(&temp);
1142
1143
/* Calculate it. */
1144
bc_divide (num1, num2, &temp, scale);
1145
if (quot)
1146
quotient = bc_copy_num (temp);
1147
bc_multiply (temp, num2, &temp, rscale);
1148
bc_sub (num1, temp, rem, rscale);
1149
bc_free_num (&temp);
1150
1151
if (quot)
1152
{
1153
bc_free_num (quot);
1154
*quot = quotient;
1155
}
1156
1157
return 0; /* Everything is OK. */
1158
}
1159
1160
1161
/* Modulo for numbers. This computes NUM1 % NUM2 and puts the
1162
result in RESULT. */
1163
1164
int
1165
bc_modulo (num1, num2, result, scale)
1166
bc_num num1, num2, *result;
1167
int scale;
1168
{
1169
return bc_divmod (num1, num2, NULL, result, scale);
1170
}
1171
1172
/* Raise BASE to the EXPO power, reduced modulo MOD. The result is
1173
placed in RESULT. If a EXPO is not an integer,
1174
only the integer part is used. */
1175
1176
int
1177
bc_raisemod (base, expo, mod, result, scale)
1178
bc_num base, expo, mod, *result;
1179
int scale;
1180
{
1181
bc_num power, exponent, parity, temp;
1182
int rscale;
1183
1184
/* Check for correct numbers. */
1185
if (bc_is_zero(mod)) return -1;
1186
if (bc_is_neg(expo)) return -1;
1187
1188
/* Set initial values. */
1189
power = bc_copy_num (base);
1190
exponent = bc_copy_num (expo);
1191
temp = bc_copy_num (_one_);
1192
bc_init_num(&parity);
1193
1194
/* Check the base for scale digits. */
1195
if (base->n_scale != 0)
1196
bc_rt_warn ("non-zero scale in base");
1197
1198
/* Check the exponent for scale digits. */
1199
if (exponent->n_scale != 0)
1200
{
1201
bc_rt_warn ("non-zero scale in exponent");
1202
bc_divide (exponent, _one_, &exponent, 0); /*truncate */
1203
}
1204
1205
/* Check the modulus for scale digits. */
1206
if (mod->n_scale != 0)
1207
bc_rt_warn ("non-zero scale in modulus");
1208
1209
/* Do the calculation. */
1210
rscale = MAX(scale, base->n_scale);
1211
while ( !bc_is_zero(exponent) )
1212
{
1213
(void) bc_divmod (exponent, _two_, &exponent, &parity, 0);
1214
if ( !bc_is_zero(parity) )
1215
{
1216
bc_multiply (temp, power, &temp, rscale);
1217
(void) bc_modulo (temp, mod, &temp, scale);
1218
}
1219
1220
bc_multiply (power, power, &power, rscale);
1221
(void) bc_modulo (power, mod, &power, scale);
1222
}
1223
1224
/* Assign the value. */
1225
bc_free_num (&power);
1226
bc_free_num (&exponent);
1227
bc_free_num (result);
1228
*result = temp;
1229
return 0; /* Everything is OK. */
1230
}
1231
1232
/* Raise NUM1 to the NUM2 power. The result is placed in RESULT.
1233
Maximum exponent is LONG_MAX. If a NUM2 is not an integer,
1234
only the integer part is used. */
1235
1236
void
1237
bc_raise (num1, num2, result, scale)
1238
bc_num num1, num2, *result;
1239
int scale;
1240
{
1241
bc_num temp, power;
1242
long exponent;
1243
int rscale;
1244
int pwrscale;
1245
int calcscale;
1246
char neg;
1247
1248
/* Check the exponent for scale digits and convert to a long. */
1249
if (num2->n_scale != 0)
1250
bc_rt_warn ("non-zero scale in exponent");
1251
exponent = bc_num2long (num2);
1252
if (exponent == 0 && (num2->n_len > 1 || num2->n_value[0] != 0))
1253
bc_rt_error ("exponent too large in raise");
1254
1255
/* Special case if exponent is a zero. */
1256
if (exponent == 0)
1257
{
1258
bc_free_num (result);
1259
*result = bc_copy_num (_one_);
1260
return;
1261
}
1262
1263
/* Other initializations. */
1264
if (exponent < 0)
1265
{
1266
neg = TRUE;
1267
exponent = -exponent;
1268
rscale = scale;
1269
}
1270
else
1271
{
1272
neg = FALSE;
1273
rscale = MIN (num1->n_scale*exponent, MAX(scale, num1->n_scale));
1274
}
1275
1276
/* Set initial value of temp. */
1277
power = bc_copy_num (num1);
1278
pwrscale = num1->n_scale;
1279
while ((exponent & 1) == 0)
1280
{
1281
pwrscale = 2*pwrscale;
1282
bc_multiply (power, power, &power, pwrscale);
1283
exponent = exponent >> 1;
1284
}
1285
temp = bc_copy_num (power);
1286
calcscale = pwrscale;
1287
exponent = exponent >> 1;
1288
1289
/* Do the calculation. */
1290
while (exponent > 0)
1291
{
1292
pwrscale = 2*pwrscale;
1293
bc_multiply (power, power, &power, pwrscale);
1294
if ((exponent & 1) == 1) {
1295
calcscale = pwrscale + calcscale;
1296
bc_multiply (temp, power, &temp, calcscale);
1297
}
1298
exponent = exponent >> 1;
1299
}
1300
1301
/* Assign the value. */
1302
if (neg)
1303
{
1304
bc_divide (_one_, temp, result, rscale);
1305
bc_free_num (&temp);
1306
}
1307
else
1308
{
1309
bc_free_num (result);
1310
*result = temp;
1311
if ((*result)->n_scale > rscale)
1312
(*result)->n_scale = rscale;
1313
}
1314
bc_free_num (&power);
1315
}
1316
1317
/* Take the square root NUM and return it in NUM with SCALE digits
1318
after the decimal place. */
1319
1320
int
1321
bc_sqrt (num, scale)
1322
bc_num *num;
1323
int scale;
1324
{
1325
int rscale, cmp_res, done;
1326
int cscale;
1327
bc_num guess, guess1, point5, diff;
1328
1329
/* Initial checks. */
1330
cmp_res = bc_compare (*num, _zero_);
1331
if (cmp_res < 0)
1332
return 0; /* error */
1333
else
1334
{
1335
if (cmp_res == 0)
1336
{
1337
bc_free_num (num);
1338
*num = bc_copy_num (_zero_);
1339
return 1;
1340
}
1341
}
1342
cmp_res = bc_compare (*num, _one_);
1343
if (cmp_res == 0)
1344
{
1345
bc_free_num (num);
1346
*num = bc_copy_num (_one_);
1347
return 1;
1348
}
1349
1350
/* Initialize the variables. */
1351
rscale = MAX (scale, (*num)->n_scale);
1352
bc_init_num(&guess);
1353
bc_init_num(&guess1);
1354
bc_init_num(&diff);
1355
point5 = bc_new_num (1,1);
1356
point5->n_value[1] = 5;
1357
1358
1359
/* Calculate the initial guess. */
1360
if (cmp_res < 0)
1361
{
1362
/* The number is between 0 and 1. Guess should start at 1. */
1363
guess = bc_copy_num (_one_);
1364
cscale = (*num)->n_scale;
1365
}
1366
else
1367
{
1368
/* The number is greater than 1. Guess should start at 10^(exp/2). */
1369
bc_int2num (&guess,10);
1370
1371
bc_int2num (&guess1,(*num)->n_len);
1372
bc_multiply (guess1, point5, &guess1, 0);
1373
guess1->n_scale = 0;
1374
bc_raise (guess, guess1, &guess, 0);
1375
bc_free_num (&guess1);
1376
cscale = 3;
1377
}
1378
1379
/* Find the square root using Newton's algorithm. */
1380
done = FALSE;
1381
while (!done)
1382
{
1383
bc_free_num (&guess1);
1384
guess1 = bc_copy_num (guess);
1385
bc_divide (*num, guess, &guess, cscale);
1386
bc_add (guess, guess1, &guess, 0);
1387
bc_multiply (guess, point5, &guess, cscale);
1388
bc_sub (guess, guess1, &diff, cscale+1);
1389
if (bc_is_near_zero (diff, cscale))
1390
{
1391
if (cscale < rscale+1)
1392
cscale = MIN (cscale*3, rscale+1);
1393
else
1394
done = TRUE;
1395
}
1396
}
1397
1398
/* Assign the number and clean up. */
1399
bc_free_num (num);
1400
bc_divide (guess,_one_,num,rscale);
1401
bc_free_num (&guess);
1402
bc_free_num (&guess1);
1403
bc_free_num (&point5);
1404
bc_free_num (&diff);
1405
return 1;
1406
}
1407
1408
1409
/* The following routines provide output for bcd numbers package
1410
using the rules of POSIX bc for output. */
1411
1412
/* This structure is used for saving digits in the conversion process. */
1413
typedef struct stk_rec {
1414
long digit;
1415
struct stk_rec *next;
1416
} stk_rec;
1417
1418
/* The reference string for digits. */
1419
static char ref_str[] = "0123456789ABCDEF";
1420
1421
1422
/* A special output routine for "multi-character digits." Exactly
1423
SIZE characters must be output for the value VAL. If SPACE is
1424
non-zero, we must output one space before the number. OUT_CHAR
1425
is the actual routine for writing the characters. */
1426
1427
void
1428
bc_out_long (val, size, space, out_char)
1429
long val;
1430
int size, space;
1431
#ifdef __STDC__
1432
void (*out_char)(int);
1433
#else
1434
void (*out_char)();
1435
#endif
1436
{
1437
char digits[40];
1438
int len, ix;
1439
1440
if (space) (*out_char) (' ');
1441
sprintf (digits, "%ld", val);
1442
len = strlen (digits);
1443
while (size > len)
1444
{
1445
(*out_char) ('0');
1446
size--;
1447
}
1448
for (ix=0; ix < len; ix++)
1449
(*out_char) (digits[ix]);
1450
}
1451
1452
/* Output of a bcd number. NUM is written in base O_BASE using OUT_CHAR
1453
as the routine to do the actual output of the characters. */
1454
1455
void
1456
bc_out_num (num, o_base, out_char, leading_zero)
1457
bc_num num;
1458
int o_base;
1459
#ifdef __STDC__
1460
void (*out_char)(int);
1461
#else
1462
void (*out_char)();
1463
#endif
1464
int leading_zero;
1465
{
1466
char *nptr;
1467
int index, fdigit, pre_space;
1468
stk_rec *digits, *temp;
1469
bc_num int_part, frac_part, base, cur_dig, t_num, max_o_digit;
1470
1471
/* The negative sign if needed. */
1472
if (num->n_sign == MINUS) (*out_char) ('-');
1473
1474
/* Output the number. */
1475
if (bc_is_zero (num))
1476
(*out_char) ('0');
1477
else
1478
if (o_base == 10)
1479
{
1480
/* The number is in base 10, do it the fast way. */
1481
nptr = num->n_value;
1482
if (num->n_len > 1 || *nptr != 0)
1483
for (index=num->n_len; index>0; index--)
1484
(*out_char) (BCD_CHAR(*nptr++));
1485
else
1486
nptr++;
1487
1488
if (leading_zero && bc_is_zero (num))
1489
(*out_char) ('0');
1490
1491
/* Now the fraction. */
1492
if (num->n_scale > 0)
1493
{
1494
(*out_char) ('.');
1495
for (index=0; index<num->n_scale; index++)
1496
(*out_char) (BCD_CHAR(*nptr++));
1497
}
1498
}
1499
else
1500
{
1501
/* special case ... */
1502
if (leading_zero && bc_is_zero (num))
1503
(*out_char) ('0');
1504
1505
/* The number is some other base. */
1506
digits = NULL;
1507
bc_init_num (&int_part);
1508
bc_divide (num, _one_, &int_part, 0);
1509
bc_init_num (&frac_part);
1510
bc_init_num (&cur_dig);
1511
bc_init_num (&base);
1512
bc_sub (num, int_part, &frac_part, 0);
1513
/* Make the INT_PART and FRAC_PART positive. */
1514
int_part->n_sign = PLUS;
1515
frac_part->n_sign = PLUS;
1516
bc_int2num (&base, o_base);
1517
bc_init_num (&max_o_digit);
1518
bc_int2num (&max_o_digit, o_base-1);
1519
1520
1521
/* Get the digits of the integer part and push them on a stack. */
1522
while (!bc_is_zero (int_part))
1523
{
1524
bc_modulo (int_part, base, &cur_dig, 0);
1525
temp = (stk_rec *) malloc (sizeof(stk_rec));
1526
if (temp == NULL) bc_out_of_memory();
1527
temp->digit = bc_num2long (cur_dig);
1528
temp->next = digits;
1529
digits = temp;
1530
bc_divide (int_part, base, &int_part, 0);
1531
}
1532
1533
/* Print the digits on the stack. */
1534
if (digits != NULL)
1535
{
1536
/* Output the digits. */
1537
while (digits != NULL)
1538
{
1539
temp = digits;
1540
digits = digits->next;
1541
if (o_base <= 16)
1542
(*out_char) (ref_str[ (int) temp->digit]);
1543
else
1544
bc_out_long (temp->digit, max_o_digit->n_len, 1, out_char);
1545
free (temp);
1546
}
1547
}
1548
1549
/* Get and print the digits of the fraction part. */
1550
if (num->n_scale > 0)
1551
{
1552
(*out_char) ('.');
1553
pre_space = 0;
1554
t_num = bc_copy_num (_one_);
1555
while (t_num->n_len <= num->n_scale) {
1556
bc_multiply (frac_part, base, &frac_part, num->n_scale);
1557
fdigit = bc_num2long (frac_part);
1558
bc_int2num (&int_part, fdigit);
1559
bc_sub (frac_part, int_part, &frac_part, 0);
1560
if (o_base <= 16)
1561
(*out_char) (ref_str[fdigit]);
1562
else {
1563
bc_out_long (fdigit, max_o_digit->n_len, pre_space, out_char);
1564
pre_space = 1;
1565
}
1566
bc_multiply (t_num, base, &t_num, 0);
1567
}
1568
bc_free_num (&t_num);
1569
}
1570
1571
/* Clean up. */
1572
bc_free_num (&int_part);
1573
bc_free_num (&frac_part);
1574
bc_free_num (&base);
1575
bc_free_num (&cur_dig);
1576
bc_free_num (&max_o_digit);
1577
}
1578
}
1579
/* Convert a number NUM to a long. The function returns only the integer
1580
part of the number. For numbers that are too large to represent as
1581
a long, this function returns a zero. This can be detected by checking
1582
the NUM for zero after having a zero returned. */
1583
1584
long
1585
bc_num2long (num)
1586
bc_num num;
1587
{
1588
long val;
1589
char *nptr;
1590
int index;
1591
1592
/* Extract the int value, ignore the fraction. */
1593
val = 0;
1594
nptr = num->n_value;
1595
for (index=num->n_len; (index>0) && (val<=(LONG_MAX/BASE)); index--)
1596
val = val*BASE + *nptr++;
1597
1598
/* Check for overflow. If overflow, return zero. */
1599
if (index>0) val = 0;
1600
if (val < 0) val = 0;
1601
1602
/* Return the value. */
1603
if (num->n_sign == PLUS)
1604
return (val);
1605
else
1606
return (-val);
1607
}
1608
1609
1610
/* Convert an integer VAL to a bc number NUM. */
1611
1612
void
1613
bc_int2num (num, val)
1614
bc_num *num;
1615
int val;
1616
{
1617
char buffer[30];
1618
char *bptr, *vptr;
1619
int ix = 1;
1620
char neg = 0;
1621
1622
/* Sign. */
1623
if (val < 0)
1624
{
1625
neg = 1;
1626
val = -val;
1627
}
1628
1629
/* Get things going. */
1630
bptr = buffer;
1631
*bptr++ = val % BASE;
1632
val = val / BASE;
1633
1634
/* Extract remaining digits. */
1635
while (val != 0)
1636
{
1637
*bptr++ = val % BASE;
1638
val = val / BASE;
1639
ix++; /* Count the digits. */
1640
}
1641
1642
/* Make the number. */
1643
bc_free_num (num);
1644
*num = bc_new_num (ix, 0);
1645
if (neg) (*num)->n_sign = MINUS;
1646
1647
/* Assign the digits. */
1648
vptr = (*num)->n_value;
1649
while (ix-- > 0)
1650
*vptr++ = *--bptr;
1651
}
1652
1653
/* Convert a numbers to a string. Base 10 only.*/
1654
1655
char
1656
*num2str (num)
1657
bc_num num;
1658
{
1659
char *str, *sptr;
1660
char *nptr;
1661
int index, signch;
1662
1663
/* Allocate the string memory. */
1664
signch = ( num->n_sign == PLUS ? 0 : 1 ); /* Number of sign chars. */
1665
if (num->n_scale > 0)
1666
str = (char *) malloc (num->n_len + num->n_scale + 2 + signch);
1667
else
1668
str = (char *) malloc (num->n_len + 1 + signch);
1669
if (str == NULL) bc_out_of_memory();
1670
1671
/* The negative sign if needed. */
1672
sptr = str;
1673
if (signch) *sptr++ = '-';
1674
1675
/* Load the whole number. */
1676
nptr = num->n_value;
1677
for (index=num->n_len; index>0; index--)
1678
*sptr++ = BCD_CHAR(*nptr++);
1679
1680
/* Now the fraction. */
1681
if (num->n_scale > 0)
1682
{
1683
*sptr++ = '.';
1684
for (index=0; index<num->n_scale; index++)
1685
*sptr++ = BCD_CHAR(*nptr++);
1686
}
1687
1688
/* Terminate the string and return it! */
1689
*sptr = '\0';
1690
return (str);
1691
}
1692
/* Convert strings to bc numbers. Base 10 only.*/
1693
1694
void
1695
bc_str2num (num, str, scale)
1696
bc_num *num;
1697
char *str;
1698
int scale;
1699
{
1700
int digits, strscale;
1701
char *ptr, *nptr;
1702
char zero_int;
1703
1704
/* Prepare num. */
1705
bc_free_num (num);
1706
1707
/* Check for valid number and count digits. */
1708
ptr = str;
1709
digits = 0;
1710
strscale = 0;
1711
zero_int = FALSE;
1712
if ( (*ptr == '+') || (*ptr == '-')) ptr++; /* Sign */
1713
while (*ptr == '0') ptr++; /* Skip leading zeros. */
1714
while (isdigit((int)*ptr)) ptr++, digits++; /* digits */
1715
if (*ptr == '.') ptr++; /* decimal point */
1716
while (isdigit((int)*ptr)) ptr++, strscale++; /* digits */
1717
if ((*ptr != '\0') || (digits+strscale == 0))
1718
{
1719
*num = bc_copy_num (_zero_);
1720
return;
1721
}
1722
1723
/* Adjust numbers and allocate storage and initialize fields. */
1724
strscale = MIN(strscale, scale);
1725
if (digits == 0)
1726
{
1727
zero_int = TRUE;
1728
digits = 1;
1729
}
1730
*num = bc_new_num (digits, strscale);
1731
1732
/* Build the whole number. */
1733
ptr = str;
1734
if (*ptr == '-')
1735
{
1736
(*num)->n_sign = MINUS;
1737
ptr++;
1738
}
1739
else
1740
{
1741
(*num)->n_sign = PLUS;
1742
if (*ptr == '+') ptr++;
1743
}
1744
while (*ptr == '0') ptr++; /* Skip leading zeros. */
1745
nptr = (*num)->n_value;
1746
if (zero_int)
1747
{
1748
*nptr++ = 0;
1749
digits = 0;
1750
}
1751
for (;digits > 0; digits--)
1752
*nptr++ = CH_VAL(*ptr++);
1753
1754
1755
/* Build the fractional part. */
1756
if (strscale > 0)
1757
{
1758
ptr++; /* skip the decimal point! */
1759
for (;strscale > 0; strscale--)
1760
*nptr++ = CH_VAL(*ptr++);
1761
}
1762
}
1763
1764
/* pn prints the number NUM in base 10. */
1765
1766
static void
1767
out_char (int c)
1768
{
1769
putchar(c);
1770
}
1771
1772
1773
void
1774
pn (num)
1775
bc_num num;
1776
{
1777
bc_out_num (num, 10, out_char, 0);
1778
out_char ('\n');
1779
}
1780
1781
1782
/* pv prints a character array as if it was a string of bcd digits. */
1783
void
1784
pv (name, num, len)
1785
char *name;
1786
unsigned char *num;
1787
int len;
1788
{
1789
int i;
1790
printf ("%s=", name);
1791
for (i=0; i<len; i++) printf ("%c",BCD_CHAR(num[i]));
1792
printf ("\n");
1793
}
Loggerhead is a web-based interface for Breezy
Version: 2.0.1