|
|
1.1 root 1: /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
2:
3: /*
4: $Header: b1num.c,v 1.4 85/08/22 16:51:59 timo Exp $
5: */
6:
7: /* B numbers, basic external interface */
8:
9: #include "b.h"
10: #include "b0con.h"
11: #include "b1obj.h"
12: #include "b1num.h"
13:
14: /*
15: * This file contains operations on numbers that are not predefined
16: * B functions (the latter are defined in `bfun.c').
17: * This includes conversion of numeric B values to C `int's and
18: * `double's, (numval() and intval()),
19: * but also utilities for comparing numeric values and hashing a numeric
20: * value to something usable as initialization for the random generator
21: * without chance of overflow (so numval(v) is unusable).
22: * It is also possible to build numbers of all types using mk_integer,
23: * mk_exact (or mk_rat) and mk_approx. Note the rather irregular
24: * (historical!) argument structure for these: mk_approx has a
25: * `double' argument, where mk_exact and mk_rat have two values
26: * which must be `integer' (not `int's!) as arguments.
27: * The random number generator used by the DRAW and CHOOSE statements
28: * is also in this file.
29: */
30:
31: /*
32: * ival is used internally by intval() and large().
33: * It converts an integer to double precision, yielding
34: * a huge value when overflow occurs (but giving no error).
35: */
36:
37: Hidden double ival(u) integer u; {
38: double x = 0;
39: register int i;
40:
41: if (IsSmallInt(u)) return SmallIntVal(u);
42: for (i = Length(u)-1; i >= 0; --i) {
43: if (x >= Maxreal/BASE)
44: return Msd(u) < 0 ? -Maxreal : Maxreal;
45: x = x*BASE + Digit(u, i);
46: }
47:
48: return x;
49: }
50:
51:
52: /* Determine if a value may be converted to an int */
53:
54: Visible bool large(v) value v; {
55: double r;
56: if (!Is_number(v) || !integral(v)) {
57: error(MESS(1300, "number not an integer"));
58: return No;
59: }
60: if (Rational(v)) v= (value) Numerator((rational)v);
61: r= ival((integer)v);
62: if (r > Maxint) return Yes;
63: return No;
64: }
65:
66: /* return the C `int' value of a B numeric value, if it exists. */
67:
68: Visible int intval(v) value v; {
69: /* v must be an Integral number or a Rational with Denominator==1
70: which may result from n round x [via mk_exact]!. */
71: double i;
72: if (IsSmallInt(v)) i= SmallIntVal(v);
73: else {
74: if (!Is_number(v)) syserr(MESS(1301, "intval on non-number"));
75: if (!integral(v)) {
76: error(MESS(1302, "number not an integer"));
77: return 0;
78: }
79: if (Rational(v)) v= (value) Numerator((rational)v);
80: i= ival((integer)v);
81: }
82: if (i > Maxint || i < -Maxint) {
83: error(MESS(1303, "exceedingly large integer"));
84: return 0;
85: }
86: return (int) i;
87: }
88:
89:
90: /* convert an int to an intlet */
91:
92: Visible intlet propintlet(i) int i; {
93: if (i > Maxintlet || i < -Maxintlet) {
94: error(MESS(1304, "exceedingly large integer"));
95: return 0;
96: }
97: return i;
98: }
99:
100:
101: /*
102: * determine if a number is integer
103: */
104:
105: Visible bool integral(v) value v; {
106: if (Integral(v) || (Rational(v) && Denominator((rational)v) == int_1))
107: return Yes;
108: else return No;
109: }
110:
111: /*
112: * mk_exact makes an exact number out of two integers.
113: * The third argument is the desired number of digits after the decimal
114: * point when the number is printed (see comments in `bfun.c' for
115: * `round' and code in `bnuC.c').
116: * This printing size (if positive) is hidden in the otherwise nearly
117: * unused * 'Size' field of the value struct
118: * (cf. definition of Rational(v) etc.).
119: */
120:
121: Visible value mk_exact(p, q, len) integer p, q; int len; {
122: rational r = mk_rat(p, q, len);
123:
124: if (Denominator(r) == int_1 && len <= 0) {
125: integer i = (integer) Copy(Numerator(r));
126: release((value) r);
127: return (value) i;
128: }
129:
130: return (value) r;
131: }
132:
133: #define uSMALL 1
134: #define uINT 2
135: #define uRAT 4
136: #define uAPP 8
137: #define vSMALL 16
138: #define vINT 32
139: #define vRAT 64
140: #define vAPP 128
141:
142: Visible relation numcomp(u, v) value u, v; {
143: int tu, tv; relation s;
144:
145: if (IsSmallInt(u)) tu = uSMALL;
146: else if (Length(u) >= 0) tu = uINT;
147: else if (Length(u) <= -2) tu = uRAT;
148: else tu = uAPP;
149: if (IsSmallInt(v)) tv = vSMALL;
150: else if (Length(v) >= 0) tv = vINT;
151: else if (Length(v) <= -2) tv = vRAT;
152: else tv = vAPP;
153:
154: switch (tu|tv) {
155:
156: case uSMALL|vSMALL: return SmallIntVal(u) - SmallIntVal(v);
157:
158: case uSMALL|vINT:
159: case uINT|vSMALL:
160: case uINT|vINT: return int_comp((integer)u, (integer)v);
161:
162: case uSMALL|vRAT:
163: case uINT|vRAT:
164: u = (value) mk_rat((integer)u, int_1, 0);
165: s = rat_comp((rational)u, (rational)v);
166: release(u);
167: return s;
168:
169: case uRAT|vRAT:
170: return rat_comp((rational)u, (rational)v);
171:
172: case uRAT|vSMALL:
173: case uRAT|vINT:
174: v = (value) mk_rat((integer)v, int_1, 0);
175: s = rat_comp((rational)u, (rational)v);
176: release(v);
177: return s;
178:
179: case uSMALL|vAPP:
180: case uINT|vAPP:
181: case uRAT|vAPP:
182: u = approximate(u);
183: s = app_comp((real)u, (real)v);
184: release(u);
185: return s == 0 ? -1 : s; /* u < ~u = v */
186:
187: case uAPP|vAPP:
188: return app_comp((real)u, (real)v);
189:
190: case uAPP|vSMALL:
191: case uAPP|vINT:
192: case uAPP|vRAT:
193: v = approximate(v);
194: s = app_comp((real)u, (real)v);
195: release(v);
196: return s == 0 ? 1 : s; /* u = ~v > v */
197:
198: default: syserr(MESS(1305, "num_comp")); /* NOTREACHED */
199:
200: }
201: }
202:
203:
204: /*
205: * Deliver 10**n, where n is a (maybe negative!) C integer.
206: * The result is a value (integer or rational, actually).
207: */
208:
209: Visible value tento(n) int n; {
210: if (n < 0) {
211: integer i= int_tento(-n);
212: value v= (value) mk_exact(int_1, i, 0);
213: release((value) i);
214: return v;
215: }
216: return (value) int_tento(n);
217: }
218:
219:
220: /*
221: * numval returns the numerical value of any numeric B value
222: * as a C `double'.
223: */
224:
225: Visible double numval(u) value u; {
226: double expo, frac;
227:
228: if (!Is_number(u)) {
229: error(MESS(1306, "value not a number"));
230: return 0.0;
231: }
232: u = approximate(u);
233: expo = Expo((real)u), frac = Frac((real)u);
234: release(u);
235: if (expo > Maxexpo) {
236: error(MESS(1307, "approximate number too large to be handled"));
237: return 0.0;
238: }
239: if(expo < Minexpo) return 0.0;
240: return ldexp(frac, (int)expo);
241: }
242:
243:
244: /*
245: * Random numbers
246: */
247:
248:
249: /*
250: * numhash produces a `double' number that depends on the value of
251: * v, useful for initializing the random generator.
252: * Needs rewriting, so that for instance numhash(n) never equals n.
253: * The magic numbers here are chosen at random.
254: */
255:
256: Visible double numhash(v) value v; {
257: if (Integral(v)) {
258: double d = 0;
259: register int i;
260:
261: if (IsSmallInt(v)) return SmallIntVal(v);
262: for (i = Length(v) - 1; i >= 0; --i) {
263: d *= 2;
264: d += Digit((integer)v, i);
265: }
266:
267: return d;
268: }
269:
270: if (Rational(v))
271: return .777 * numhash((value) Numerator((rational)v)) +
272: .123 * numhash((value) Denominator((rational)v));
273:
274: return numval(v); /* Which fails for HUGE reals. Never mind. */
275: }
276:
277:
278: /* Initialize the random generator */
279:
280: double lastran;
281:
282: Hidden Procedure setran (seed) double seed; {
283: double x;
284:
285: x = seed >= 0 ? seed : -seed;
286: while (x >= 1) x /= 10;
287: lastran = floor(67108864.0*x);
288: }
289:
290: Visible Procedure set_random(v) value v; {
291: setran((double) hash(v));
292: }
293:
294:
295: /* Return a random number in [0, 1). */
296:
297: Visible value random() {
298: double p;
299:
300: p = 26353589.0 * lastran + 1;
301: lastran = p - 67108864.0*floor(p/67108864.0);
302:
303: return (value) mk_approx(lastran / 67108864.0, 0.0);
304: }
305:
306: Visible Procedure initnum() {
307: rat_init();
308: setran((double) SEED);
309: }
310:
311: Visible Procedure endnum() {
312: endrat();
313: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.