|
|
1.1 root 1: /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1984. */
2: /* $Header: B1num.c,v 1.1 84/06/28 00:48:56 timo Exp $ */
3:
4: /* B numbers, small version */
5:
6: /*
7: * THIS VERSION SHOULD ONLY BE USED IF
8: * THE SYSTEM IS TOO LARGE OTHERWISE.
9: * IT USES FLOATING POINT ARITHMETIC FOR EXACT NUMBERS
10: * INSTEAD OF ARBITRARY LENGTH RATIONAL ARITHMETIC.
11: */
12:
13: #include "b.h"
14: #include "b0con.h"
15: #include "b1obj.h"
16: #include "b2syn.h" /* for def of Keymark() only */
17: #include "B1num.h"
18:
19: value numerator(v) value v; {
20: Checknum(v);
21: if (!Exact(v)) error("*/ on approximate number");
22: return mk_int(Numerator(v));
23: }
24:
25: value denominator(v) value v; {
26: Checknum(v);
27: if (!Exact(v)) error("/* on approximate number"); /* */
28: return mk_int(Denominator(v));
29: }
30:
31: double numval(v) value v; {
32: Checknum(v);
33: return Numval(v);
34: }
35:
36: checkint(v) value v; {
37: Checknum(v);
38: if (Denominator(v) != One) error("number not an integer");
39: }
40:
41: bool large(v) value v; {
42: checkint(v);
43: if (Numerator(v) < -Maxint || Numerator(v) > Maxint) return Yes;
44: return No;
45: }
46:
47: int intval(v) value v; {
48: checkint(v);
49: return (int)Numerator(v);
50: }
51:
52: intlet propintlet(i) int i; {
53: if (i < -Maxintlet || i > Maxintlet)
54: error("exceedingly large integer");
55: return i;
56: }
57:
58: integer gcd(i, j) integer i, j; {
59: integer k;
60: if (i == Zero && j == Zero) syserr("gcd(0, 0)");
61: if (i != floor(i) || j != floor(j))
62: syserr("gcd called with non-integer");
63: if (i < Zero) i= -i; if (j < Zero) j= -j;
64: if (i < j) {
65: k= i; i= j; j= k;
66: }
67: while (j >= One) {
68: k= i-j*floor(i/j);
69: i= j; j= k;
70: }
71: if (j != Zero) error(
72: "arithmetic overflow while simplifying exact number");
73: if (i != floor(i)) syserr("gcd returns non-integer");
74: return i;
75: }
76:
77: value b_zero, b_one, b_minus_one, zero, one;
78:
79: value mk_exact(p, q, len) register integer p, q; intlet len; {
80: value v; integer d;
81: if (q == One && len ==0) {
82: if (p == Zero) return copy(b_zero);
83: if (p == One) return copy(b_one);
84: if (p == -One) return copy(b_minus_one);
85: }
86: v= grab_num(len);
87: if (q == One) {
88: Numerator(v)= p; Denominator(v)= q;
89: return v;
90: }
91: if (q == Zero) error("attempt to make exact number with denominator 0");
92: if (q < Zero) {p= -p; q= -q;}
93: d= (q == One ? One : p == One ? One : gcd(p, q));
94: Numerator(v)= p/d; Denominator(v)= q/d;
95: return v;
96: }
97:
98: bool integral(v) value v; {
99: return Integral(v);
100: }
101:
102: value mk_integer(p) int p; {
103: return mk_exact((integer)p, One, 0);
104: }
105:
106: value mk_int(p) integer p; {
107: return mk_exact(p, One, 0);
108: }
109:
110: value mk_approx(x) register double x; {
111: value v= grab_num(0);
112: Approxval(v)= x; Denominator(v)= Zero;
113: return v;
114: }
115:
116: initnum() {
117: b_zero= grab_num(0);
118: Numerator(b_zero)= Zero; Denominator(b_zero)= One;
119: b_one= grab_num(0);
120: Numerator(b_one)= One; Denominator(b_one)= One;
121: b_minus_one= grab_num(0);
122: Numerator(b_minus_one)= -One; Denominator(b_minus_one)= One;
123: zero= mk_integer(0);
124: one= mk_integer(1);
125: }
126:
127: value approximate(v) value v; {
128: if (!Exact(v)) return copy(v);
129: return mk_approx(Numerator(v)/Denominator(v));
130: }
131:
132: numcomp(v, w) value v, w; {
133: double vv= Numval(v), ww= Numval(w);
134: if (vv < ww) return -1;
135: if (vv > ww) return 1;
136: if (Exact(v) && Exact(w)) return 0;
137: if (Exact(v)) return -1; /* 1 < 1E0 */
138: if (Exact(w)) return 1; /* 1E0 > 1 */
139: return 0;
140: }
141:
142: double numhash(v) value v; {
143: number *n= (number *)Ats(v);
144: return .123*n->p + .777*n->q;
145: }
146:
147: #define CONVBUFSIZ 100
148: char convbuf[CONVBUFSIZ];
149:
150: string convnum(v) value v; {
151: double x; string bp; bool prec_loss= No;
152: Checknum(v);
153: x= Numval(v);
154: conv: if (!prec_loss && Exact(v) && fabs(x) <= LONG &&
155: fabs(Numerator(v)) < BIG && fabs(Denominator(v)) < BIG) {
156: intlet len= 0 < Length(v) && Length(v) <= MAXNUMDIG ? Length(v) : 0;
157: intlet dcnt, sigcnt; bool sig;
158: if (Denominator(v) != One) {
159: intlet k; double p= 1.0, q;
160: prec_loss= Yes;
161: for (k= 1; k < MAXNUMDIG; k++) {
162: p*= 10.0;
163: q= p/Denominator(v);
164: if (k >= len && q == floor(q)) {
165: prec_loss= No;
166: break;
167: }
168: }
169: len= k;
170: }
171: convex: sprintf(convbuf, "%.*f", len, x);
172: dcnt= sigcnt= 0; sig= No;
173: for (bp= convbuf; *bp != '\0'; bp++)
174: if ('0' <= *bp && *bp <= '9') {
175: dcnt++;
176: if (*bp != '0') sig= Yes;
177: if (sig) sigcnt++;
178: }
179: if (sigcnt < MINNUMDIG && prec_loss) goto conv;
180: if (dcnt > MAXNUMDIG) {
181: if (len <= 0) syserr("conversion error 1");
182: if (Denominator(v) == One) len= 0;
183: else len-= dcnt-MAXNUMDIG;
184: if (len < 0) syserr("conversion error 2");
185: goto convex;
186: }
187: } else { /*approx etc*/
188: sprintf(convbuf, "%.*e", MAXNUMDIG-5, x);
189: for (bp= convbuf; *bp != '\0'; bp++)
190: if (*bp == 'e') {
191: *bp= 'E';
192: break;
193: }
194: }
195: return convbuf;
196: }
197:
198: value numconst(tx, q) txptr tx, q; {
199: bool dig= No; double ex= 0, ap= 1; intlet ndap, len= 0;
200: while (tx < q && '0' <= *tx && *tx <= '9') {
201: dig= Yes;
202: ex= 10*ex+(*tx++ - '0');
203: }
204: if (tx < q && *tx == '.') {
205: tx++; ndap= 0;
206: while (tx < q && '0' <= *tx && *tx <= '9') {
207: dig= Yes; ndap++;
208: len= *tx == '0' ? ndap : 0;
209: ex= 10*ex+(*tx++ - '0'); ap*= 10;
210: }
211: if (!dig) syserr("numconst[1]");
212: }
213: if (tx < q && *tx == 'E') {
214: intlet sign= 1; double expo= 0;
215: tx++;
216: if (!('0' <= *tx && *tx <= '9') && Keymark(*tx)) {
217: tx--;
218: goto exact;
219: }
220: if (!dig) ex= 1;
221: if (tx < q && (*tx == '+' || *tx == '-'))
222: if (*tx++ == '-') sign= -1;
223: dig= No;
224: while (tx < q && '0' <= *tx && *tx <= '9') {
225: dig= Yes;
226: expo= 10*expo+(*tx++ - '0');
227: }
228: if (!dig) syserr("numconst[2]");
229: return mk_approx(ex/ap*exp(sign*expo*log(10.0)));
230: }
231: exact: return mk_exact(ex, ap, len);
232: }
233:
234: printnum(f1, v) FILE *f1; value v; {
235: FILE *f= f1 ? f1 : stdout;
236: if (!Exact(v) || Denominator(v) == One) {
237: if (!Exact(v))
238: fputc('~', f);
239: fputs(convnum(v), f);
240: }
241: else {
242: value w = numerator(v);
243: fputs(convnum(w), f);
244: release(w);
245: fputc('/', f);
246: w = denominator(v);
247: fputs(convnum(w), f);
248: release(w);
249: }
250: if (!f1) fputc('\n', f); /* Flush buffer for sdb */
251: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.