1 /****************************************************************/
2 /* file arith.c
3
4 ARIBAS interpreter for Arithmetic
5 Copyright (C) 1996 O.Forster
6
7 This program is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2 of the License, or
10 (at your option) any later version.
11
12 This program is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with this program; if not, write to the Free Software
19 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
20
21 Address of the author
22
23 Otto Forster
24 Math. Institut der LMU
25 Theresienstr. 39
26 D-80333 Muenchen, Germany
27
28 Email forster@mathematik.uni-muenchen.de
29 */
30 /****************************************************************/
31 /*
32 ** arith.c
33 ** arithmetic functions
34 **
35 ** date of last change
36 ** 1995-03-11: set_floatprec, get_floatprec
37 ** 1997-02-20: fixed bug in Frandom()
38 ** 1997-02-25: modified Fmod()
39 ** 1997-03-31: error message in Gcompare
40 ** 1997-04-13: reorg (newintsym)
41 ** 2001-01-06: multiprec floats
42 ** 2002-02-24: vector arithmetic, 1st step
43 ** 2002-04-19: divide(x,y)
44 ** 2003-02-11: addition of gf2n_int's
45 ** 2003-02-15: Fdecode made safe wrt garbage collection
46 ** 2003-11-21: removed bug in vecdiv
47 ** 2007-02-23: changed negatevec to make it robust wrt gc
48 ** 2007-11-15: function roundhalf: case of number exactly = 1/2
49 */
50
51 #include "common.h"
52
53 PUBLIC void iniarith _((void));
54 PUBLIC truc addints _((truc *ptr, int minflg));
55 PUBLIC truc scalintvec _((truc *ptr1, truc *ptr2));
56 PUBLIC unsigned random2 _((unsigned u));
57 PUBLIC unsigned random4 _((unsigned u));
58 PUBLIC int cmpnums _((truc *ptr1, truc *ptr2, int type));
59
60
61 PUBLIC truc zero, constone;
62
63 PUBLIC truc sfloatsym, dfloatsym, lfloatsym, xfloatsym;
64 PUBLIC truc realsym, integsym, int_sym;
65 PUBLIC truc plussym, minussym, uminsym;
66 PUBLIC truc timessym, divsym, modsym, divfsym, powersym;
67 PUBLIC truc ariltsym, arigtsym, arilesym, arigesym;
68
69 PUBLIC long maxfltex, maxdecex, exprange;
70
71 /*-------------------------------------------------------------*/
72
73 PRIVATE truc Fplus _((void));
74 PRIVATE truc addfloats _((truc *ptr, int minflg));
75 PRIVATE truc Fminus _((void));
76 PRIVATE truc Fnegate _((void));
77 PRIVATE truc Sinc _((void));
78 PRIVATE truc Sdec _((void));
79 PRIVATE truc Sincaux _((truc symb, int s));
80 PRIVATE truc Fabsolute _((void));
81 PRIVATE truc Fmax _((int argn));
82 PRIVATE truc Fmin _((int argn));
83 PRIVATE truc Fmaxint _((void));
84 PRIVATE truc Gminmax _((truc symb, int argn));
85 PRIVATE truc Fodd _((void));
86 PRIVATE truc Feven _((void));
87 PRIVATE int odd _((truc *ptr));
88 PRIVATE truc Ftimes _((void));
89 PRIVATE truc multints _((truc *ptr));
90 PRIVATE truc multfloats _((truc *ptr));
91 PRIVATE truc Fsum _((void));
92 PRIVATE truc sumintvec _((truc *argptr, int argn));
93 PRIVATE truc sumfltvec _((truc *argptr, int argn));
94 PRIVATE truc Fprod _((void));
95 PRIVATE truc prodintvec _((truc *argptr, int argn));
96 PRIVATE truc prodfloatvec _((truc *argptr, int argn));
97 PRIVATE truc Fdivf _((void));
98 PRIVATE truc Fdiv _((void));
99 PRIVATE truc Fdivide _((void));
100 PRIVATE truc Fmod _((void));
101 PRIVATE truc Gvecmod _((int flg));
102 PRIVATE truc divide _((truc *ptr, int tflag));
103 PRIVATE truc modout _((truc *ptr));
104 PRIVATE truc divfloats _((truc *ptr));
105 PRIVATE truc Ftrunc _((void));
106 PRIVATE truc Fround _((void));
107 PRIVATE truc Ffloor _((void));
108 PRIVATE truc Ffrac _((void));
109 PRIVATE truc Gtruncaux _((truc symb));
110 PRIVATE void intfrac _((numdata *npt1, numdata *np2));
111 PRIVATE int roundhalf _((numdata *nptr));
112 PRIVATE void floshiftint _((numdata *nptr));
113 PRIVATE truc Fpower _((void));
114 PRIVATE truc exptints _((truc *ptr, unsigned a));
115 PRIVATE truc exptfloats _((truc *ptr));
116 PRIVATE int exptovfl _((word2 *x, int n, unsigned a));
117 PRIVATE truc Fisqrt _((void));
118 PRIVATE int cmpfloats _((truc *ptr1, truc *ptr2, int prec));
119 PRIVATE truc Farilt _((void));
120 PRIVATE truc Farigt _((void));
121 PRIVATE truc Farile _((void));
122 PRIVATE truc Farige _((void));
123 PRIVATE int Gcompare _((truc symb));
124 PRIVATE void inirandstate _((word2 *rr));
125 PRIVATE void nextrand _((word2 *rr, int n));
126 PRIVATE truc Frandom _((void));
127 PRIVATE truc Frandseed _((int argn));
128 PRIVATE int objfltprec _((truc obj));
129 PRIVATE truc Ffloat _((int argn));
130 PRIVATE truc Fsetfltprec _((int argn));
131 PRIVATE truc Fgetfltprec _((int argn));
132 PRIVATE truc Fmaxfltprec _((void));
133 PRIVATE truc Fsetprnprec _((void));
134 PRIVATE truc Fgetprnprec _((void));
135 PRIVATE int precdesc _((truc obj));
136 PRIVATE truc Fdecode _((void));
137 PRIVATE truc Fbitnot _((void));
138 PRIVATE truc Fbitset _((void));
139 PRIVATE truc Fbitclear _((void));
140 PRIVATE truc Gbitset _((truc symb));
141 PRIVATE truc Fbittest _((void));
142 PRIVATE truc Fbitshift _((void));
143 PRIVATE truc Fbitlength _((void));
144 PRIVATE truc Fbitcount _((void));
145 PRIVATE truc Fcardinal _((void));
146 PRIVATE truc Finteger _((void));
147 PRIVATE truc Gcardinal _((truc symb));
148 PRIVATE truc Fbitand _((void));
149 PRIVATE truc Fbitor _((void));
150 PRIVATE truc Fbitxor _((void));
151 PRIVATE truc Gboole _((truc symb, ifunaa boolfun));
152
153 PRIVATE truc negatevec _((truc *ptr));
154 PRIVATE truc addvecs _((truc sym, truc *ptr));
155 PRIVATE truc addintvecs _((truc *ptr1, truc *ptr2));
156 PRIVATE truc addfltvecs _((truc *ptr1, truc *ptr2));
157 PRIVATE truc scalvec _((truc *ptr1, truc *ptr2));
158 PRIVATE truc scalfltvec _((truc *ptr1, truc *ptr2));
159 PRIVATE truc vecdiv _((truc *vptr, truc *zz));
160 PRIVATE truc vecdivfloat _((truc *vptr, truc *zz));
161
162 PRIVATE int chkplusargs _((truc sym, truc *argptr));
163 PRIVATE int chktimesargs _((truc *argptr));
164 PRIVATE int chkmodargs _((truc sym, truc *argptr));
165 PRIVATE int chkdivfargs _((truc sym, truc *argptr));
166
167 PRIVATE truc floatsym;
168 PRIVATE truc decodsym;
169 PRIVATE truc sumsym, prodsym;
170 PRIVATE truc bnotsym, bandsym, borsym, bxorsym, bshiftsym;
171 PRIVATE truc btestsym, bsetsym, bclrsym, blensym, bcountsym;
172 PRIVATE truc cardsym;
173
174 PRIVATE truc maxsym, minsym;
175 PRIVATE truc maxintsym;
176 PRIVATE truc setfltsym, getfltsym, maxfltsym, setprnsym, getprnsym;
177 PRIVATE truc incsym, decsym;
178 PRIVATE truc abssym, oddsym, evensym;
179 PRIVATE truc isqrtsym;
180 PRIVATE truc dividesym, div_sym, mod_sym;
181 PRIVATE truc truncsym, roundsym, fracsym, floorsym;
182 PRIVATE truc randsym, rseedsym;
183
184 PRIVATE word2 RandSeed[4];
185 PRIVATE word4 MaxBits;
186 PRIVATE int curFltPrec;
187
188 /*--------------------------------------------------------*/
iniarith()189 PUBLIC void iniarith()
190 {
191 word4 u;
192
193 zero = mkfixnum(0);
194 constone = mkfixnum(1);
195
196 integsym = newsym("integer", sTYPESPEC, zero);
197 int_sym = new0symsig("integer", sFBINARY, (wtruc)Finteger, s_ii);
198 realsym = newsym("real", sTYPESPEC, fltzero(deffltprec()));
199
200 sfloatsym = newsym("single_float", sSYMBCONST, mkfixnum(FltPrec[0]));
201 dfloatsym = newsym("double_float", sSYMBCONST, mkfixnum(FltPrec[1]));
202 lfloatsym = newsym("long_float", sSYMBCONST, mkfixnum(FltPrec[2]));
203 xfloatsym = newsym("extended_float",sSYMBCONST, mkfixnum(FltPrec[3]));
204
205 floatsym = newsymsig("float", sFBINARY, (wtruc)Ffloat, s_12rn);
206 setfltsym = newsymsig("set_floatprec", sFBINARY,
207 (wtruc)Fsetfltprec, s_12);
208 getfltsym = newsymsig("get_floatprec", sFBINARY,
209 (wtruc)Fgetfltprec, s_01);
210 maxfltsym = newsymsig("max_floatprec", sFBINARY, (wtruc)Fmaxfltprec,s_0);
211 setprnsym = newsymsig("set_printprec", sFBINARY, (wtruc)Fsetprnprec,s_1);
212 getprnsym = newsymsig("get_printprec", sFBINARY, (wtruc)Fgetprnprec,s_0);
213
214 decodsym = newsymsig("decode_float",sFBINARY,(wtruc)Fdecode, s_vr);
215 plussym = newintsym("+", sFBINARY, (wtruc)Fplus);
216 minussym = newintsym("-", sFBINARY, (wtruc)Fminus);
217 uminsym = newintsym("-", sFBINARY, (wtruc)Fnegate);
218 timessym = newintsym("*", sFBINARY, (wtruc)Ftimes);
219 divfsym = newintsym("/", sFBINARY, (wtruc)Fdivf);
220 powersym = newintsym("**", sFBINARY, (wtruc)Fpower);
221 sumsym = newsymsig("sum", sFBINARY, (wtruc)Fsum, s_nv);
222 prodsym = newsymsig("product",sFBINARY, (wtruc)Fprod, s_nv);
223
224 divsym = newintsym("div",sFBINARY, (wtruc)Fdiv);
225 div_sym = newsym("div", sINFIX, divsym);
226 SYMcc1(div_sym) = DIVTOK;
227 dividesym = newsymsig("divide",sFBINARY, (wtruc)Fdivide, s_2);
228
229 modsym = newintsym("mod",sFBINARY, (wtruc)Fmod);
230 mod_sym = newsym("mod", sINFIX, modsym);
231 SYMcc1(mod_sym) = MODTOK;
232
233 abssym = newsymsig("abs", sFBINARY, (wtruc)Fabsolute,s_rr);
234 maxsym = newsymsig("max", sFBINARY, (wtruc)Fmax,s_1u);
235 minsym = newsymsig("min", sFBINARY, (wtruc)Fmin,s_1u);
236 maxintsym = newsymsig("max_intsize", sFBINARY, (wtruc)Fmaxint,s_0);
237
238 oddsym = newsymsig("odd", sFBINARY, (wtruc)Fodd,s_ii);
239 evensym = newsymsig("even", sFBINARY, (wtruc)Feven,s_ii);
240 incsym = newsymsig("inc", sSBINARY, (wtruc)Sinc,s_12ii);
241 decsym = newsymsig("dec", sSBINARY, (wtruc)Sdec,s_12ii);
242 isqrtsym = newsymsig("isqrt", sFBINARY, (wtruc)Fisqrt,s_ii);
243
244 truncsym = newsymsig("trunc", sFBINARY, (wtruc)Ftrunc, s_rr);
245 roundsym = newsymsig("round", sFBINARY, (wtruc)Fround, s_rr);
246 floorsym = newsymsig("floor", sFBINARY, (wtruc)Ffloor, s_rr);
247 fracsym = newsymsig("frac", sFBINARY, (wtruc)Ffrac, s_rr);
248
249 randsym = newsymsig("random", sFBINARY, (wtruc)Frandom, s_rr);
250 rseedsym = newsymsig("random_seed",sFBINARY,(wtruc)Frandseed, s_01);
251
252 ariltsym = newintsym("<", sFBINARY, (wtruc)Farilt);
253 arigtsym = newintsym(">", sFBINARY, (wtruc)Farigt);
254 arilesym = newintsym("<=", sFBINARY, (wtruc)Farile);
255 arigesym = newintsym(">=", sFBINARY, (wtruc)Farige);
256
257 bnotsym = newsymsig("bit_not", sFBINARY, (wtruc)Fbitnot, s_ii);
258 bandsym = newsymsig("bit_and", sFBINARY, (wtruc)Fbitand, s_iii);
259 borsym = newsymsig("bit_or", sFBINARY, (wtruc)Fbitor, s_iii);
260 bxorsym = newsymsig("bit_xor", sFBINARY, (wtruc)Fbitxor, s_iii);
261 btestsym = newsymsig("bit_test", sFBINARY, (wtruc)Fbittest,s_iii);
262 bsetsym = newsymsig("bit_set", sFBINARY, (wtruc)Fbitset, s_iii);
263 bclrsym = newsymsig("bit_clear",sFBINARY, (wtruc)Fbitclear, s_iii);
264 bshiftsym = newsymsig("bit_shift",sFBINARY, (wtruc)Fbitshift, s_iii);
265 blensym = newsymsig("bit_length",sFBINARY,(wtruc)Fbitlength,s_ii);
266 bcountsym = newsymsig("bit_count",sFBINARY, (wtruc)Fbitcount,s_ii);
267 cardsym = newsymsig("cardinal", sFBINARY, (wtruc)Fcardinal,s_ii);
268
269 u = aribufSize;
270 u <<= 4;
271 MaxBits = u;
272 u -= 256;
273 if(u > MAXFLTLIM)
274 u = MAXFLTLIM;
275 maxfltex = u;
276 maxdecex = (u/10) * 3;
277 exprange = u - u/3 + u/38; /* log(2) = 1 - 1/3 + 1/38 */
278
279 inirandstate(RandSeed);
280 iniaritx();
281 iniarity();
282 iniaritz();
283 }
284 /*--------------------------------------------------------*/
285 #define DIVFLAG 1
286 #define MODFLAG 2
287 #define DDIVFLAG (DIVFLAG | MODFLAG)
288 /*--------------------------------------------------------*/
Fplus()289 PRIVATE truc Fplus()
290 {
291 truc res;
292 int type;
293
294 type = chkplusargs(plussym,argStkPtr-1);
295 if(type > fBIGNUM) {
296 curFltPrec = deffltprec();
297 res = addfloats(argStkPtr-1,0);
298 }
299 else switch(type) {
300 case fFIXNUM:
301 case fBIGNUM:
302 res = addints(argStkPtr-1,0);
303 break;
304 case fVECTOR:
305 res = addvecs(plussym,argStkPtr-1);
306 break;
307 case fGF2NINT:
308 res = addgf2ns(argStkPtr-1);
309 break;
310 case aERROR:
311 default:
312 res = brkerr();
313 break;
314 }
315 return(res);
316 }
317 /*--------------------------------------------------------*/
addints(ptr,minflag)318 PUBLIC truc addints(ptr,minflag)
319 truc *ptr;
320 int minflag; /* if minflag != 0, subtract */
321 {
322 word2 *y;
323 int n1, n2, n;
324 int sign1, sign2, sign, cmp;
325
326 n1 = bigretr(ptr,AriBuf,&sign1);
327 n2 = bigref(ptr+1,&y,&sign2);
328 if(minflag) {
329 sign2 = (sign2 ? 0 : MINUSBYTE);
330 }
331 if(sign1 == sign2) {
332 sign = sign1;
333 n = addarr(AriBuf,n1,y,n2);
334 }
335 else {
336 cmp = cmparr(AriBuf,n1,y,n2);
337 if(cmp > 0) {
338 sign = sign1;
339 n = subarr(AriBuf,n1,y,n2);
340 }
341 else if(cmp < 0) {
342 sign = sign2;
343 n = sub1arr(AriBuf,n1,y,n2);
344 }
345 else
346 return(zero);
347 }
348 return(mkint(sign,AriBuf,n));
349 }
350 /*-------------------------------------------------------------*/
addfloats(ptr,minflag)351 PRIVATE truc addfloats(ptr,minflag)
352 truc *ptr;
353 int minflag;
354 {
355 numdata accum, temp;
356 int prec, cmp;
357
358 prec = curFltPrec + 1;
359 if(prec < 3)
360 prec++;
361 accum.digits = AriBuf;
362 temp.digits = AriScratch;
363 getnumalign(prec,ptr,&accum);
364 getnumalign(prec,ptr+1,&temp);
365 if(minflag)
366 temp.sign = (temp.sign ? 0 : MINUSBYTE);
367 adjustoffs(&accum,&temp);
368 if(accum.sign == temp.sign) {
369 accum.len =
370 addarr(accum.digits,accum.len,temp.digits,temp.len);
371 }
372 else {
373 cmp = cmparr(accum.digits,accum.len,temp.digits,temp.len);
374 if(cmp > 0) {
375 accum.len =
376 subarr(accum.digits,accum.len,temp.digits,temp.len);
377 }
378 else if(cmp < 0) {
379 accum.sign = temp.sign;
380 accum.len =
381 sub1arr(accum.digits,accum.len,temp.digits,temp.len);
382 }
383 else {
384 accum.len = 0;
385 }
386 }
387 return(mkfloat(curFltPrec,&accum));
388 }
389 /*--------------------------------------------------------*/
Fminus()390 PRIVATE truc Fminus()
391 {
392 truc res;
393 int type;
394
395 type = chkplusargs(minussym,argStkPtr-1);
396 if(type > fBIGNUM) {
397 curFltPrec = deffltprec();
398 res = addfloats(argStkPtr-1,-1);
399 }
400 else switch(type) {
401 case fFIXNUM:
402 case fBIGNUM:
403 res = addints(argStkPtr-1,-1);
404 break;
405 case fVECTOR:
406 res = addvecs(minussym,argStkPtr-1);
407 break;
408 case fGF2NINT:
409 res = addgf2ns(argStkPtr-1);
410 break;
411 case aERROR:
412 default:
413 res = brkerr();
414 break;
415 }
416 return(res);
417 }
418 /*--------------------------------------------------------*/
Fnegate()419 PRIVATE truc Fnegate()
420 {
421 int flg;
422 truc res[1];
423
424 flg = *FLAGPTR(argStkPtr);
425
426 if(flg >= fFIXNUM) {
427 res[0] = mkcopy(argStkPtr);
428 changesign(res);
429 return res[0];
430 }
431 else if(flg == fVECTOR) {
432 flg = chknumvec(uminsym,argStkPtr);
433 if(flg == aERROR)
434 return brkerr();
435 else
436 return negatevec(argStkPtr);
437 }
438 else if(flg == fGF2NINT) {
439 return *argStkPtr;
440 }
441 else { /* flg == aERROR */
442 error(uminsym,err_num,*argStkPtr);
443 return brkerr();
444 }
445 }
446 /*--------------------------------------------------------*/
Sinc()447 PRIVATE truc Sinc()
448 {
449 return(Sincaux(incsym,1));
450 }
451 /*--------------------------------------------------------*/
Sdec()452 PRIVATE truc Sdec()
453 {
454 return(Sincaux(decsym,-1));
455 }
456 /*--------------------------------------------------------*/
Sincaux(symb,s)457 PRIVATE truc Sincaux(symb,s)
458 truc symb;
459 int s;
460 {
461 truc res;
462 long number;
463 int argn;
464 int flg;
465
466 argn = *ARGCOUNTPTR(evalStkPtr);
467 res = eval(ARG1PTR(evalStkPtr));
468 ARGpush(res);
469 if(argn == 2) {
470 res = eval(ARGNPTR(evalStkPtr,2));
471 ARGpush(res);
472 }
473 flg = chkints(symb,argStkPtr-argn+1,argn);
474 if(flg == aERROR) {
475 ARGnpop(argn);
476 return(brkerr());
477 }
478 else if(argn == 1 && flg == fFIXNUM) {
479 number = *WORD2PTR(argStkPtr);
480 if(*SIGNPTR(argStkPtr))
481 number = -number;
482 res = mkinum(number + s);
483 ARGpop();
484 }
485 else {
486 if(argn == 1)
487 ARGpush(constone);
488 s = (s > 0 ? 0 : -1);
489 res = addints(argStkPtr-1,s);
490 ARGnpop(2);
491 }
492 return(Lvalassign(ARG1PTR(evalStkPtr),res));
493 }
494 /*--------------------------------------------------------*/
Ftimes()495 PRIVATE truc Ftimes()
496 {
497 truc res;
498 int type;
499
500 type = chktimesargs(argStkPtr-1);
501 if(type <= fBIGNUM) {
502 switch(type) {
503 case fFIXNUM:
504 case fBIGNUM:
505 res = multints(argStkPtr-1);
506 break;
507 case fGF2NINT:
508 res = multgf2ns(argStkPtr-1);
509 break;
510 default: /* type == aERROR */
511 res = brkerr();
512 }
513 }
514 else if((type & 0xFF00) == 0) { /* float obj */
515 curFltPrec = deffltprec();
516 res = multfloats(argStkPtr-1);
517 }
518 else if((type >> 8) == fVECTOR) {
519 type |= 0xFF;
520 if(type >= fFIXNUM) {
521 res = scalvec(argStkPtr-1,argStkPtr);
522 }
523 else if(type == fGF2NINT) {
524 error(timessym,err_imp,voidsym);
525 res = brkerr();
526 }
527 else
528 res = brkerr();
529 }
530 else {
531 error(timessym,err_case,mksfixnum(type));
532 res = brkerr();
533 }
534 return(res);
535 }
536 /*----------------------------------------------------------*/
537 /*
538 ** multiply integers in ptr[0] and ptr[1]
539 */
multints(ptr)540 PRIVATE truc multints(ptr)
541 truc *ptr;
542 {
543 word2 *x, *y;
544 int n1, n2, n, sign, sign2;
545
546 n1 = bigref(ptr,&x,&sign);
547 n2 = bigref(ptr+1,&y,&sign2);
548 if(n1 + n2 >= aribufSize)
549 goto errexit;
550 else if(!n1 || !n2)
551 return(zero);
552 n = multbig(x,n1,y,n2,AriBuf,AriScratch);
553 sign = (sign == sign2 ? 0 : MINUSBYTE);
554 return(mkint(sign,AriBuf,n));
555 errexit:
556 error(timessym,err_ovfl,voidsym);
557 return(brkerr());
558 }
559 /*---------------------------------------------------------*/
multfloats(ptr)560 PRIVATE truc multfloats(ptr)
561 truc *ptr;
562 {
563 numdata prod, temp;
564 int prec;
565 int n;
566
567 prec = curFltPrec + 1;
568 prod.digits = AriBuf;
569 n = getnumtrunc(prec,ptr,&prod);
570 refnumtrunc(prec,ptr+1,&temp);
571 n = multtrunc(prec,&prod,&temp,AriScratch);
572
573 if(n < 0) {
574 error(timessym,err_ovfl,voidsym);
575 return(brkerr());
576 }
577 return(mkfloat(curFltPrec,&prod));
578 }
579 /*--------------------------------------------------------*/
negatevec(ptr)580 PRIVATE truc negatevec(ptr)
581 truc *ptr;
582 {
583 truc res;
584 truc tmp[1];
585 int len,k;
586
587 WORKpush(zero);
588 WORKpush(mkcopy(ptr));
589 len = *VECLENPTR(workStkPtr);
590 for(k=0; k<len; k++) {
591 workStkPtr[-1] = *(VECTORPTR(workStkPtr)+k);
592 tmp[0] = mkcopy(workStkPtr - 1);
593 *(VECTORPTR(workStkPtr)+k) = changesign(tmp);
594 }
595 res = WORKretr();
596 WORKpop();
597 return res;
598 }
599 /*--------------------------------------------------------*/
600 /*
601 ** addition of two vectors (int or float) ptr[0] and ptr[1]
602 */
addvecs(sym,ptr)603 PRIVATE truc addvecs(sym,ptr)
604 truc sym;
605 truc *ptr;
606 {
607 truc tmp;
608 int len1, len2, flg1, flg;
609
610 flg1 = chknumvec(sym,ptr);
611 if(flg1 == aERROR)
612 return brkerr();
613 flg = chknumvec(sym,ptr+1);
614 if(flg == aERROR) {
615 return brkerr();
616 }
617 if(flg < flg1)
618 flg = flg1;
619
620 if(sym == minussym) {
621 ptr[1] = mkcopy(ptr+1);
622 ptr[1] = negatevec(ptr+1);
623 }
624 len1 = *VECLENPTR(ptr);
625 len2 = *VECLENPTR(ptr+1);
626 if(len2 < len1) {
627 tmp = mkcopy(ptr);
628 ptr[0] = ptr[1];
629 ptr[1] = tmp;
630 }
631 else if(sym == plussym) {
632 ptr[1] = mkcopy(ptr+1);
633 }
634 if(flg <= fBIGNUM) {
635 return addintvecs(ptr,ptr+1);
636 }
637 else {
638 curFltPrec = deffltprec();
639 return addfltvecs(ptr,ptr+1);
640 }
641 }
642 /*--------------------------------------------------------*/
643 /*
644 ** adds two integer vectors *ptr1 and *ptr2
645 ** length of *ptr1 must be <= length of *ptr2
646 */
addintvecs(ptr1,ptr2)647 PRIVATE truc addintvecs(ptr1,ptr2)
648 truc *ptr1,*ptr2;
649 {
650 truc obj;
651 int len, k;
652
653 len = *VECLENPTR(ptr1);
654 WORKnpush(2);
655 for(k=0; k<len; k++) {
656 workStkPtr[-1] = *(VECTORPTR(ptr1) + k);
657 workStkPtr[0] = *(VECTORPTR(ptr2) + k);
658 obj = addints(workStkPtr-1,0);
659 *(VECTORPTR(ptr2) + k) = obj;
660 }
661 WORKnpop(2);
662 return *ptr2;
663 }
664 /*--------------------------------------------------------*/
665 /*
666 ** adds two float vectors *ptr1 and *ptr2
667 ** length of *ptr1 must be <= length of *ptr2
668 */
addfltvecs(ptr1,ptr2)669 PRIVATE truc addfltvecs(ptr1,ptr2)
670 truc *ptr1,*ptr2;
671 {
672 truc obj;
673 int len, k;
674
675 len = *VECLENPTR(ptr1);
676 WORKnpush(2);
677 for(k=0; k<len; k++) {
678 workStkPtr[-1] = *(VECTORPTR(ptr1) + k);
679 workStkPtr[0] = *(VECTORPTR(ptr2) + k);
680 obj = addfloats(workStkPtr-1,0);
681 *(VECTORPTR(ptr2) + k) = obj;
682 }
683 WORKnpop(2);
684 return *ptr2;
685 }
686 /*---------------------------------------------------------*/
687 /*
688 ** multiplies vector given in *ptr2 by scalar given in *ptr1
689 */
scalvec(ptr1,ptr2)690 PRIVATE truc scalvec(ptr1,ptr2)
691 truc *ptr1, *ptr2;
692 {
693 int flg, flg0;
694
695 flg = chknumvec(timessym,ptr2);
696 if(flg == aERROR)
697 return brkerr();
698 flg0 = *FLAGPTR(ptr1);
699 if(flg0 > flg)
700 flg = flg0;
701 if(flg <= fBIGNUM) {
702 return scalintvec(ptr1,ptr2);
703 }
704 else {
705 curFltPrec = deffltprec();
706 return scalfltvec(ptr1,ptr2);
707 }
708 }
709 /*---------------------------------------------------------*/
710 /*
711 ** Multiplies int vector given in *ptr2 by scalar in *ptr1
712 ** uses AriBuf and AriScratch
713 */
scalintvec(ptr1,ptr2)714 PUBLIC truc scalintvec(ptr1,ptr2)
715 truc *ptr1, *ptr2;
716 {
717 truc obj;
718 int len, k;
719
720 WORKpush(*ptr1);
721 WORKpush(constone);
722 obj = mkcopy(ptr2);
723 WORKpush(obj);
724 len = *VECLENPTR(ptr2);
725 for(k=0; k<len; k++) {
726 *(workStkPtr-1) = *(VECTORPTR(workStkPtr) + k);
727 obj = multints(workStkPtr-2);
728 *(VECTORPTR(workStkPtr) + k) = obj;
729 }
730 obj = WORKretr();
731 WORKnpop(2);
732 return obj;
733 }
734 /*---------------------------------------------------------*/
735 /*
736 ** Multiplies float vector given in *ptr2 by scalar in *ptr1
737 */
scalfltvec(ptr1,ptr2)738 PRIVATE truc scalfltvec(ptr1,ptr2)
739 truc *ptr1, *ptr2;
740 {
741 truc obj;
742 int len, k;
743
744 WORKpush(*ptr1);
745 WORKpush(constone);
746 obj = mkcopy(ptr2);
747 WORKpush(obj);
748 len = *VECLENPTR(ptr2);
749 for(k=0; k<len; k++) {
750 *(workStkPtr-1) = *(VECTORPTR(workStkPtr) + k);
751 obj = multfloats(workStkPtr-2);
752 *(VECTORPTR(workStkPtr) + k) = obj;
753 }
754 obj = WORKretr();
755 WORKnpop(2);
756 return obj;
757 }
758 /*--------------------------------------------------------*/
Fabsolute()759 PRIVATE truc Fabsolute()
760 {
761 if(chknum(abssym,argStkPtr) == aERROR)
762 return(brkerr());
763
764 *argStkPtr = mkcopy(argStkPtr);
765 wipesign(argStkPtr);
766 return(*argStkPtr);
767 }
768 /*--------------------------------------------------------*/
Fmax(argn)769 PRIVATE truc Fmax(argn)
770 int argn;
771 {
772 return(Gminmax(maxsym,argn));
773 }
774 /*--------------------------------------------------------*/
Fmin(argn)775 PRIVATE truc Fmin(argn)
776 int argn;
777 {
778 return(Gminmax(minsym,argn));
779 }
780 /*--------------------------------------------------------*/
Fmaxint()781 PRIVATE truc Fmaxint()
782 {
783 long mxint = maxdecex - 8;
784
785 return(mkinum(mxint));
786 }
787 /*--------------------------------------------------------*/
Gminmax(symb,argn)788 PRIVATE truc Gminmax(symb,argn)
789 truc symb;
790 int argn;
791 {
792 struct vector *vptr;
793 truc *ptr, *argptr;
794 int cmp, type;
795
796 if(argn == 1 && *FLAGPTR(argStkPtr) == fVECTOR) {
797 vptr = (struct vector *)TAddress(argStkPtr);
798 argn = vptr->len;
799 if(argn == 0) {
800 error(symb,err_args,*argStkPtr);
801 goto errexit;
802 }
803 argptr = &(vptr->ele0);
804 }
805 else
806 argptr = argStkPtr - argn + 1;
807 if((type = chknums(symb,argptr,argn)) == aERROR)
808 goto errexit;
809 ptr = argptr++;
810 while(--argn > 0) {
811 cmp = cmpnums(ptr,argptr,type);
812 if(symb == minsym)
813 cmp = -cmp;
814 if(cmp < 0)
815 ptr = argptr;
816 argptr++;
817 }
818 return(*ptr);
819 errexit:
820 return brkerr();
821 }
822 /*--------------------------------------------------------*/
Fodd()823 PRIVATE truc Fodd()
824 {
825 int ret = odd(argStkPtr);
826
827 if(ret == aERROR)
828 return(brkerr());
829 return(ret ? true : false);
830 }
831 /*--------------------------------------------------------*/
Feven()832 PRIVATE truc Feven()
833 {
834 int ret = odd(argStkPtr);
835
836 if(ret == aERROR)
837 return(brkerr());
838 return(ret ? false : true);
839 }
840 /*--------------------------------------------------------*/
odd(ptr)841 PRIVATE int odd(ptr)
842 truc *ptr;
843 {
844 word2 *x;
845 int sign;
846
847 if(bigref(ptr,&x,&sign) == aERROR)
848 return(aERROR);
849 else
850 return(x[0] & 1);
851 }
852 /*---------------------------------------------------------*/
Fsum()853 PRIVATE truc Fsum()
854 {
855 struct vector *vec;
856 truc *ptr;
857 long res;
858 int flg;
859 int len;
860
861 if(*FLAGPTR(argStkPtr) != fVECTOR) {
862 error(sumsym,err_vect,*argStkPtr);
863 return(brkerr());
864 }
865 vec = (struct vector *)TAddress(argStkPtr);
866 len = vec->len;
867 ptr = &(vec->ele0);
868 flg = chknums(sumsym,ptr,len);
869
870 if(flg == fFIXNUM) {
871 res = 0;
872 while(--len >= 0) {
873 if(*SIGNPTR(ptr))
874 res -= *WORD2PTR(ptr);
875 else
876 res += *WORD2PTR(ptr);
877 ptr++;
878 }
879 return(mkinum(res));
880 }
881 else if(flg == fBIGNUM) {
882 return(sumintvec(ptr,len));
883 }
884 else if(flg >= fFLTOBJ) {
885 curFltPrec = deffltprec();
886 return(sumfltvec(ptr,len));
887 }
888 else /* vector elements are not numbers */
889 return(brkerr());
890 }
891 /*-------------------------------------------------------------*/
892 /*
893 ** sum up all elements of integer vector *argptr of length len
894 */
sumintvec(argptr,len)895 PRIVATE truc sumintvec(argptr,len)
896 truc *argptr;
897 int len;
898 {
899 word2 *y;
900 int sign, cmp, n0, n1, m;
901
902 n0 = n1 = 0;
903 while(--len >= 0) {
904 m = bigref(argptr,&y,&sign);
905 if(sign)
906 n1 = addarr(AriScratch,n1,y,m);
907 else
908 n0 = addarr(AriBuf,n0,y,m);
909 argptr++;
910 }
911 cmp = cmparr(AriBuf,n0,AriScratch,n1);
912 if(cmp < 0) {
913 sign = MINUSBYTE;
914 m = sub1arr(AriBuf,n0,AriScratch,n1);
915 }
916 else {
917 sign = 0;
918 m = subarr(AriBuf,n0,AriScratch,n1);
919 }
920 return(mkint(sign,AriBuf,m));
921 }
922 /*-------------------------------------------------------------*/
923 /*
924 ** sum up all elements of float vector *argptr of length len
925 */
sumfltvec(argptr,len)926 PRIVATE truc sumfltvec(argptr,len)
927 truc *argptr;
928 int len;
929 {
930 numdata accum, negsum, temp;
931 numdata *nptr;
932 int prec;
933 int cmp, sign;
934
935 prec = curFltPrec + 1;
936 accum.digits = AriBuf;
937 negsum.digits = AriScratch;
938 temp.digits = AriScratch + aribufSize;
939 accum.len = 0;
940 accum.expo = MOSTNEGEX;
941 negsum.len = 0;
942 negsum.expo = MOSTNEGEX;
943 while(--len >= 0) {
944 if(!getnumalign(prec,argptr++,&temp))
945 continue;
946 sign = temp.sign;
947 nptr = (sign ? &negsum : &accum);
948 adjustoffs(nptr,&temp);
949 nptr->len =
950 addarr(nptr->digits,nptr->len,temp.digits,temp.len);
951 }
952 adjustoffs(&accum,&negsum);
953 cmp = cmparr(accum.digits,accum.len,negsum.digits,negsum.len);
954 if(cmp < 0) {
955 accum.sign = MINUSBYTE;
956 accum.len =
957 sub1arr(accum.digits,accum.len,negsum.digits,negsum.len);
958 }
959 else {
960 accum.sign = 0;
961 accum.len =
962 subarr(accum.digits,accum.len,negsum.digits,negsum.len);
963 }
964 return(mkfloat(curFltPrec,&accum));
965 }
966 /*---------------------------------------------------------*/
Fprod()967 PRIVATE truc Fprod()
968 {
969 struct vector *vec;
970 truc *ptr;
971 int flg;
972 int len;
973
974 if(*FLAGPTR(argStkPtr) != fVECTOR) {
975 error(prodsym,err_vect,*argStkPtr);
976 return(brkerr());
977 }
978 vec = (struct vector *)TAddress(argStkPtr);
979 len = vec->len;
980
981 if(!len)
982 return(constone);
983 ptr = &(vec->ele0);
984 flg = chknums(prodsym,ptr,len);
985 if(flg == aERROR)
986 return(brkerr());
987 if(flg <= fBIGNUM)
988 return(prodintvec(ptr,len));
989 else { /* flg >= fFLTOBJ */
990 curFltPrec = deffltprec();
991 return(prodfloatvec(ptr,len));
992 }
993 }
994 /*----------------------------------------------------------*/
995 /*
996 ** multiplies all elements of integer vector *argptr of length len
997 */
prodintvec(argptr,len)998 PRIVATE truc prodintvec(argptr,len)
999 truc *argptr;
1000 int len;
1001 {
1002 word2 *y, *hilf;
1003 int n1, n, sign, sign1;
1004 unsigned a;
1005
1006 if(len == 0)
1007 return(constone);
1008
1009 n = bigref(argptr++,&y,&sign);
1010 cpyarr(y,n,AriBuf);
1011 hilf = AriScratch + aribufSize;
1012
1013 while(--len > 0) {
1014 if(*FLAGPTR(argptr) == fFIXNUM) {
1015 a = *WORD2PTR(argptr);
1016 n = multarr(AriBuf,n,a,AriBuf);
1017 if(n >= aribufSize)
1018 goto errexit;
1019 sign1 = *SIGNPTR(argptr);
1020 }
1021 else {
1022 n1 = bigref(argptr,&y,&sign1);
1023 if(n + n1 >= aribufSize)
1024 goto errexit;
1025 cpyarr(AriBuf,n,AriScratch);
1026 n = multbig(AriScratch,n,y,n1,AriBuf,hilf);
1027 }
1028 if(n == 0) {
1029 sign = 0;
1030 break;
1031 }
1032 if(sign1)
1033 sign = (sign ? 0 : MINUSBYTE);
1034 argptr++;
1035 }
1036 return(mkint(sign,AriBuf,n));
1037 errexit:
1038 error(prodsym,err_ovfl,voidsym);
1039 return(brkerr());
1040 }
1041 /*---------------------------------------------------------*/
1042 /*
1043 ** multiplies all elements of float vector *argptr of length len
1044 */
prodfloatvec(argptr,len)1045 PRIVATE truc prodfloatvec(argptr,len)
1046 truc *argptr;
1047 int len; /* len >= 1 */
1048 {
1049 numdata prod, temp;
1050 int prec;
1051 int n;
1052
1053 prec = curFltPrec + 1;
1054 prod.digits = AriBuf;
1055 n = getnumtrunc(prec,argptr++,&prod);
1056 while(--len > 0 && n > 0) {
1057 refnumtrunc(prec,argptr++,&temp);
1058 n = multtrunc(prec,&prod,&temp,AriScratch);
1059 }
1060 if(n < 0) {
1061 error(prodsym,err_ovfl,voidsym);
1062 return(brkerr());
1063 }
1064 return(mkfloat(curFltPrec,&prod));
1065 }
1066 /*---------------------------------------------------------*/
Fdivf()1067 PRIVATE truc Fdivf()
1068 {
1069 int type, flg;
1070
1071 type = chkdivfargs(divfsym,argStkPtr-1);
1072 if(type == aERROR) {
1073 return brkerr();
1074 }
1075 else if(type == fGF2NINT) {
1076 return divgf2ns(argStkPtr-1);
1077 }
1078 if((*argStkPtr == zero) ||
1079 ((flg = *FLAGPTR(argStkPtr)) >= fFLTOBJ && (flg & FLTZEROBIT))) {
1080 error(divfsym,err_div,voidsym);
1081 return brkerr();
1082 }
1083 curFltPrec = deffltprec();
1084 if((type & 0xFF00) == 0) {
1085 return(divfloats(argStkPtr-1));
1086 }
1087 /* else first argument is vector */
1088 flg = chknumvec(divfsym,argStkPtr-1);
1089 if(flg == aERROR) {
1090 return brkerr();
1091 }
1092 else {
1093 return(vecdivfloat(argStkPtr-1,argStkPtr));
1094 }
1095 }
1096 /*----------------------------------------------------------*/
Fdiv()1097 PRIVATE truc Fdiv()
1098 {
1099 int type, flg;
1100
1101 type = chkmodargs(divsym,argStkPtr-1);
1102 if(type == aERROR) {
1103 return(brkerr());
1104 }
1105 if(*argStkPtr == zero) {
1106 error(divsym,err_div,voidsym);
1107 return(brkerr());
1108 }
1109 if((type & 0xFF00) == 0)
1110 return(divide(argStkPtr-1,DIVFLAG));
1111
1112 /* else first argument is vector */
1113 flg = chkintvec(divsym,argStkPtr-1);
1114 if(flg == aERROR) {
1115 return(brkerr());
1116 }
1117 else {
1118 return(vecdiv(argStkPtr-1,argStkPtr));
1119 }
1120 }
1121 /*----------------------------------------------------------*/
1122 /*
1123 ** Divide integer vector given in *vptr by integer *zz
1124 */
vecdiv(vptr,zz)1125 PRIVATE truc vecdiv(vptr,zz)
1126 truc *vptr;
1127 truc *zz;
1128 {
1129 truc *ptr;
1130 truc obj;
1131 word2 *z, *hilf;
1132 int k, len, len1, len2, rlen, sign1, sign2;
1133
1134 hilf = AriScratch + aribufSize;
1135 len = *VECLENPTR(vptr);
1136 WORKpush(mkcopy(vptr));
1137 for(k=0; k<len; k++) {
1138 ptr = VECTORPTR(workStkPtr)+k;
1139 len1 = bigretr(ptr,AriScratch,&sign1);
1140 len2 = bigref(zz,&z,&sign2);
1141 len1 = divbig(AriScratch,len1,z,len2,AriBuf,&rlen,hilf);
1142 if(sign1 == sign2) {
1143 sign1 = 0;
1144 }
1145 else {
1146 sign1 = MINUSBYTE;
1147 if(rlen)
1148 len1 = incarr(AriBuf,len1,1);
1149 }
1150 obj = mkint(sign1,AriBuf,len1);
1151 *(VECTORPTR(workStkPtr)+k) = obj;
1152 }
1153 return(WORKretr());
1154 }
1155 /*----------------------------------------------------------*/
1156 /*
1157 ** Divide integer vector given in *vptr by number *zz
1158 */
vecdivfloat(vptr,zz)1159 PRIVATE truc vecdivfloat(vptr,zz)
1160 truc *vptr;
1161 truc *zz;
1162 {
1163 truc obj;
1164 int len, k;
1165
1166 WORKpush(zero);
1167 WORKpush(*zz);
1168 obj = mkcopy(vptr);
1169 WORKpush(obj);
1170 len = *VECLENPTR(vptr);
1171 for(k=0; k<len; k++) {
1172 *(workStkPtr-2) = *(VECTORPTR(workStkPtr) + k);
1173 obj = divfloats(workStkPtr-2);
1174 *(VECTORPTR(workStkPtr) + k) = obj;
1175 }
1176 obj = WORKretr();
1177 WORKnpop(2);
1178 return obj;
1179 }
1180 /*----------------------------------------------------------*/
Fmod()1181 PRIVATE truc Fmod()
1182 {
1183 int type, flg;
1184
1185 type = chkmodargs(modsym,argStkPtr-1);
1186 if(type == aERROR) {
1187 return(brkerr());
1188 }
1189 if(*argStkPtr == zero) {
1190 error(modsym,err_div,voidsym);
1191 return(brkerr());
1192 }
1193 if((type & 0xFF00) == 0) /* then both arguments are integers */
1194 return(modout(argStkPtr-1));
1195
1196 /* else first argument is vector */
1197 flg = chkintvec(modsym,argStkPtr-1);
1198 if(flg == aERROR) {
1199 return(brkerr());
1200 }
1201 type &= 0xFF;
1202 /* type == fFIXNUM or fBIGNUM */
1203 argStkPtr[-1] = mkcopy(argStkPtr-1);
1204 return Gvecmod(type);
1205 }
1206 /*----------------------------------------------------------*/
modout(ptr)1207 PRIVATE truc modout(ptr)
1208 truc *ptr;
1209 {
1210 word2 *y;
1211 int rlen, n1, n2, sign1, sign2;
1212
1213 n1 = bigretr(ptr,AriBuf,&sign1);
1214 n2 = bigref(ptr+1,&y,&sign2);
1215
1216 rlen = modbig(AriBuf,n1,y,n2,AriScratch);
1217 if(rlen && (sign1 != sign2))
1218 rlen = sub1arr(AriBuf,rlen,y,n2);
1219
1220 return(mkint(sign2,AriBuf,rlen));
1221 }
1222 /*----------------------------------------------------------*/
1223 /*
1224 ** mod out (destructively) vector given in argStkPtr[-1]
1225 ** by fixnum resp. bignum in argStkPtr[0]
1226 */
Gvecmod(flg)1227 PRIVATE truc Gvecmod(flg)
1228 int flg; /* flg == fFIXNUM or fBIGNUM */
1229 {
1230 truc *ptr;
1231 truc obj;
1232 word2 *zz, *arr;
1233 int k, len, len1, len2, rlen, sign1, sign2;
1234 unsigned z, x;
1235 int rest;
1236
1237 len = *VECLENPTR(argStkPtr-1);
1238 if(flg == fFIXNUM) {
1239 z = *WORD2PTR(argStkPtr);
1240 sign2 = *SIGNPTR(argStkPtr);
1241 ptr = VECTORPTR(argStkPtr-1);
1242 for(k=0; k<len; k++) {
1243 len1 = bigref(ptr,&arr,&sign1);
1244 x = modarr(arr,len1,z);
1245 rest = ((sign1 == sign2) || (x==0) ? x : z-x);
1246 if(sign2)
1247 rest = -rest;
1248 *ptr++ = mksfixnum(rest);
1249 }
1250 }
1251 else { /* flg == fBIGNUM */
1252 len2 = *BIGLENPTR(argStkPtr);
1253 sign2 = *SIGNUMPTR(argStkPtr);
1254 for(k=0; k<len; k++) {
1255 ptr = VECTORPTR(argStkPtr-1) + k;
1256 len1 = bigretr(ptr,AriBuf,&sign1);
1257 zz = BIGNUMPTR(argStkPtr);
1258 rlen = modbig(AriBuf,len1,zz,len2,AriScratch);
1259 if(rlen && (sign1 != sign2))
1260 rlen = sub1arr(AriBuf,rlen,zz,len2);
1261 obj = mkint(sign2,AriBuf,rlen);
1262 *(VECTORPTR(argStkPtr-1)+k) = obj;
1263 }
1264 }
1265 return argStkPtr[-1];
1266 }
1267 /*----------------------------------------------------------*/
Fdivide()1268 PRIVATE truc Fdivide()
1269 {
1270 int type;
1271
1272 type = chkmodargs(dividesym,argStkPtr-1);
1273 if(type == aERROR) {
1274 return(brkerr());
1275 }
1276 if(*argStkPtr == zero) {
1277 error(dividesym,err_div,voidsym);
1278 return(brkerr());
1279 }
1280 if((type & 0xFF00) == 0) /* then both arguments are integers */
1281 return(divide(argStkPtr-1,DDIVFLAG));
1282 else {/* else first argument is vector */
1283 error(dividesym,err_num,argStkPtr[-1]);
1284 return(brkerr());
1285 }
1286 }
1287 /*----------------------------------------------------------*/
divide(ptr,tflag)1288 PRIVATE truc divide(ptr,tflag)
1289 truc *ptr;
1290 int tflag; /* one of DIVFLAG, MODFLAG, DDIVFLAG */
1291 {
1292 truc res, vec;
1293 truc *vptr;
1294 word2 *y, *hilf;
1295 int len, rlen, n1, n2, sign1, sign2;
1296
1297 n1 = bigretr(ptr,AriScratch,&sign1);
1298 n2 = bigref(ptr+1,&y,&sign2);
1299 hilf = AriScratch + aribufSize;
1300 len = divbig(AriScratch,n1,y,n2,AriBuf,&rlen,hilf);
1301 if(tflag & DIVFLAG) {
1302 if(sign1 != sign2) {
1303 if(rlen)
1304 len = incarr(AriBuf,len,1);
1305 sign1 = MINUSBYTE;
1306 }
1307 else
1308 sign1 = 0;
1309 }
1310 if(tflag & MODFLAG) {
1311 if(rlen && (sign1 != sign2))
1312 rlen = sub1arr(AriScratch,rlen,y,n2);
1313 y = AriBuf + len;
1314 cpyarr(AriScratch,rlen,y);
1315 }
1316 if(tflag == DIVFLAG) {
1317 return mkint(sign1,AriBuf,len);
1318 }
1319 else if(tflag == MODFLAG) {
1320 return mkint(sign2,y,rlen);
1321 }
1322 else { /* tflag == DDIVFLAG */
1323 res = mkint(sign1,AriBuf,len);
1324 WORKpush(res);
1325 res = mkint(sign2,y,rlen);
1326 WORKpush(res);
1327 vec = mkvect0(2);
1328 vptr = VECTOR(vec);
1329 vptr[1] = WORKretr();
1330 vptr[0] = WORKretr();
1331 return vec;
1332 }
1333 }
1334 /*----------------------------------------------------------*/
1335 /*
1336 ** divides ptr[0] by ptr[1].
1337 ** Hypothesis: ptr[1] is not equal to 0.
1338 */
divfloats(ptr)1339 PRIVATE truc divfloats(ptr)
1340 truc *ptr;
1341 {
1342 numdata quot, temp;
1343 int prec;
1344 int n, m;
1345
1346 prec = curFltPrec + 1;
1347 quot.digits = AriBuf;
1348
1349 n = getnumtrunc(prec,ptr,");
1350 m = refnumtrunc(prec,ptr+1,&temp);
1351 n = divtrunc(prec,",&temp,AriScratch);
1352 if(n < 0) {
1353 error(divfsym,err_ovfl,voidsym);
1354 return(brkerr());
1355 }
1356 return(mkfloat(curFltPrec,"));
1357 }
1358 /*------------------------------------------------------------------*/
Ftrunc()1359 PRIVATE truc Ftrunc()
1360 {
1361 return Gtruncaux(truncsym);
1362 }
1363 /*------------------------------------------------------------------*/
Fround()1364 PRIVATE truc Fround()
1365 {
1366 return Gtruncaux(roundsym);
1367 }
1368 /*------------------------------------------------------------------*/
Ffloor()1369 PRIVATE truc Ffloor()
1370 {
1371 return Gtruncaux(floorsym);
1372 }
1373 /*------------------------------------------------------------------*/
Ffrac()1374 PRIVATE truc Ffrac()
1375 {
1376 return Gtruncaux(fracsym);
1377 }
1378 /*------------------------------------------------------------------*/
Gtruncaux(symb)1379 PRIVATE truc Gtruncaux(symb)
1380 truc symb;
1381 {
1382 numdata x1, y;
1383 int type;
1384 int prec, s1;
1385 int rhalf;
1386
1387 type = chknum(symb,argStkPtr);
1388 if(type == aERROR)
1389 return(brkerr());
1390 curFltPrec = fltprec(type);
1391 prec = curFltPrec + 1;
1392
1393 x1.digits = AriBuf;
1394 y.digits = AriScratch;
1395
1396 if(type <= fBIGNUM) {
1397 if(symb == truncsym || symb == floorsym)
1398 return(*argStkPtr);
1399 else if(symb == fracsym) {
1400 x1.len = 0;
1401 return(mkfloat(curFltPrec,&x1));
1402 }
1403 }
1404 /* argument is a float */
1405 getnumtrunc(prec,argStkPtr,&x1);
1406 s1 = x1.sign;
1407 intfrac(&x1,&y);
1408 if(symb == fracsym) {
1409 cpyarr(y.digits,y.len,AriBuf);
1410 y.digits = AriBuf;
1411 y.sign = s1;
1412 return(mkfloat(curFltPrec,&y));
1413 }
1414 floshiftint(&x1);
1415 if(symb == roundsym) {
1416 rhalf = roundhalf(&y);
1417 if(rhalf == 2) {
1418 if((x1.digits[0] & 1) && (x1.len))
1419 rhalf = 1;
1420 else
1421 rhalf = 0;
1422 }
1423 if(rhalf == 1) {
1424 x1.len = incarr(x1.digits,x1.len,1);
1425 x1.sign = s1;
1426 }
1427 }
1428 else if((symb == floorsym) && s1 && y.len) {
1429 x1.sign = s1;
1430 x1.len = incarr(x1.digits,x1.len,1);
1431 }
1432 return(mkint(x1.sign,x1.digits,x1.len));
1433 }
1434 /*------------------------------------------------------------------*/
1435 /*
1436 ** Von der in npt1 gegebenen Zahl wird destruktiv der ganzzahlige
1437 ** Anteil gebildet. Der Rest wird in npt2 gespeichert.
1438 ** Rest erhaelt dasselbe Vorzeichen wie die Ausgangszahl.
1439 */
intfrac(npt1,npt2)1440 PRIVATE void intfrac(npt1,npt2)
1441 numdata *npt1, *npt2;
1442 {
1443 word2 *x;
1444 long expo;
1445 int n, k;
1446
1447 n = npt1->len;
1448 if(n == 0 || npt1->expo >= 0) {
1449 int2numdat(0,npt2);
1450 return;
1451 }
1452 n = alignfloat(n+1,npt1);
1453 expo = npt1->expo;
1454 k = (-expo) >> 4; /* div 16, geht auf */
1455 if(k >= n) {
1456 cpynumdat(npt1,npt2);
1457 int2numdat(0,npt1);
1458 }
1459 else {
1460 x = npt1->digits;
1461 while(k > 0 && x[k-1] == 0)
1462 k--;
1463 cpyarr(x,k,npt2->digits);
1464 setarr(x,k,0);
1465 npt2->len = k;
1466 npt2->expo = expo;
1467 npt2->sign = npt1->sign;
1468 }
1469 }
1470 /*------------------------------------------------------------------*/
1471 /*
1472 ** Ergibt 1 oder 0, je nachdem die in *nptr gegebene Zahl absolut
1473 ** groesser oder kleiner 1/2 ist.
1474 ** If the number equals exactly 1/2, the return value is 2
1475 */
roundhalf(nptr)1476 PRIVATE int roundhalf(nptr)
1477 numdata *nptr;
1478 {
1479 long nn;
1480 word2 *xx;
1481 int i, bb, len;
1482
1483 len = nptr->len;
1484 if (len == 0)
1485 return 0;
1486 xx = nptr->digits;
1487 nn = bit_length(xx,len);
1488 if (nn < -nptr->expo)
1489 return 0;
1490 if (nn == -nptr->expo) {
1491 for (i=0; i < len-1; i++) {
1492 if (xx[i])
1493 return 1;
1494 }
1495 bb = (nn - 1) % 16;
1496 if (xx[len-1] != (1 << bb))
1497 return 1;
1498 else
1499 return 2;
1500 }
1501 else
1502 return 1;
1503 }
1504 /*------------------------------------------------------------------*/
floshiftint(nptr)1505 PRIVATE void floshiftint(nptr)
1506 numdata *nptr;
1507 {
1508 nptr->len = lshiftarr(nptr->digits,nptr->len,nptr->expo);
1509 nptr->expo = 0;
1510 }
1511 /*------------------------------------------------------------------*/
Fpower()1512 PRIVATE truc Fpower()
1513 {
1514 truc res;
1515 int flg, flg2, sign, m;
1516 unsigned int a;
1517 word2 *aa;
1518
1519 flg = *FLAGPTR(argStkPtr-1);
1520 flg2 = *FLAGPTR(argStkPtr);
1521 if(flg >= fFIXNUM) {
1522 if(flg2 >= fFIXNUM)
1523 flg = (flg2 >= flg ? flg2 : flg);
1524 else {
1525 error(powersym,err_num,*argStkPtr);
1526 flg = aERROR;
1527 }
1528 }
1529 else if(flg == fGF2NINT) {
1530 if(flg2 == fFIXNUM || flg2 == fBIGNUM) {
1531 return exptgf2n(argStkPtr-1);
1532 }
1533 else {
1534 error(powersym,err_int,*argStkPtr);
1535 flg = aERROR;
1536 }
1537 }
1538 else {
1539 error(powersym,err_intt,argStkPtr[-1]);
1540 flg = aERROR;
1541 }
1542 if(flg == aERROR)
1543 res = brkerr();
1544 else if(flg <= fBIGNUM) {
1545 if(argStkPtr[-1] == zero)
1546 res = (argStkPtr[0] == zero ? constone : zero);
1547 else if(argStkPtr[-1] == constone)
1548 res = constone;
1549 else {
1550 m = bigref(argStkPtr,&aa,&sign);
1551 if(sign) {
1552 flg = fFLTOBJ; /* vorlaeufig */
1553 }
1554 else if(m <= 2) {
1555 a = big2long(aa,m);
1556 res = exptints(argStkPtr-1,a);
1557 }
1558 else {
1559 error(powersym,err_ovfl,voidsym);
1560 res = brkerr();
1561 }
1562 }
1563 }
1564 if(flg >= fFLTOBJ) {
1565 curFltPrec = deffltprec();
1566 res = exptfloats(argStkPtr-1);
1567 }
1568 return(res);
1569 }
1570 /*----------------------------------------------------------------*/
exptints(ptr,a)1571 PRIVATE truc exptints(ptr,a)
1572 truc *ptr;
1573 unsigned a;
1574 {
1575 word2 *x, *temp, *hilf;
1576 int sign, n;
1577
1578 n = bigref(ptr,&x,&sign);
1579 if(exptovfl(x,n,a)) {
1580 error(powersym,err_ovfl,voidsym);
1581 return(brkerr());
1582 }
1583 temp = AriScratch;
1584 hilf = temp + aribufSize;
1585 n = power(x,n,a,AriBuf,temp,hilf);
1586 sign = (a & 1 ? sign : 0);
1587 return(mkint(sign,AriBuf,n));
1588 }
1589 /*----------------------------------------------------------------*/
exptfloats(ptr)1590 PRIVATE truc exptfloats(ptr)
1591 truc *ptr;
1592 {
1593 numdata acc, acc2;
1594 word2 *hilf;
1595 int prec, flg, len;
1596 int odd, sign = 0;
1597
1598 acc.digits = AriBuf;
1599 acc2.digits = AriScratch;
1600 hilf = AriScratch + aribufSize;
1601
1602 prec = curFltPrec + 1;
1603 len = getnumtrunc(prec,ptr,&acc);
1604 if(!len) {
1605 sign = numposneg(ptr+1);
1606 if(sign < 0) {
1607 goto errexit;
1608 }
1609 else {
1610 int2numdat((sign > 0 ? 0 : 1),&acc);
1611 goto ausgang;
1612 }
1613 }
1614 if(acc.sign) {
1615 flg = *FLAGPTR(ptr+1);
1616 if(flg > fBIGNUM) {
1617 goto errexit;
1618 }
1619 else {
1620 odd = (flg == fFIXNUM ?
1621 *WORD2PTR(ptr+1) & 1 :
1622 *BIGNUMPTR(ptr+1) & 1);
1623 sign = (odd ? MINUSBYTE : 0);
1624 }
1625 acc.sign = 0;
1626 }
1627 len = lognum(prec,&acc,hilf);
1628 if(len >= 0) {
1629 getnumtrunc(prec,ptr+1,&acc2);
1630 len = multtrunc(prec,&acc,&acc2,hilf);
1631 }
1632 if(len >= 0)
1633 len = expnum(prec,&acc,hilf);
1634 if(len == aERROR) {
1635 error(powersym,err_ovfl,voidsym);
1636 return(brkerr());
1637 }
1638 acc.sign = sign;
1639 ausgang:
1640 return(mkfloat(curFltPrec,&acc));
1641 errexit:
1642 error(powersym,err_pnum,*ptr);
1643 return(brkerr());
1644 }
1645 /*----------------------------------------------------------------*/
1646 /*
1647 ** Stellt fest, ob (x,n)**a zu overflow fuehrt.
1648 ** Rueckgabe aERROR bei overflow, oder 0 sonst
1649 */
exptovfl(x,n,a)1650 PRIVATE int exptovfl(x,n,a)
1651 word2 *x;
1652 int n;
1653 unsigned a;
1654 {
1655 numdata pow;
1656 word4 bitbound, b;
1657
1658 if(n == 0 || a <= 1)
1659 return(0);
1660 bitbound = MaxBits/a;
1661 b = n - 1;
1662 b <<= 4;
1663 b += bitlen(x[n-1]);
1664 if(b <= bitbound)
1665 return(0);
1666 else if(n > 1 || b > bitbound+1)
1667 return(aERROR);
1668 else if(n == 1 && x[0] == 2 && a < MaxBits)
1669 return(0);
1670 pow.digits = AriBuf;
1671 pwrtrunc(2,x[0],256,&pow,AriScratch);
1672 /******* oder mit log1_16() ****************/
1673 b = pow.expo;
1674 b += ((pow.len-1)<<4) + bitlen(pow.digits[pow.len-1]);
1675 bitbound = (MaxBits << 8)/a;
1676 return(b > bitbound ? aERROR : 0);
1677 }
1678 /*------------------------------------------------------------------*/
Fisqrt()1679 PRIVATE truc Fisqrt()
1680 {
1681 word2 *x, *z, *hilf;
1682 int sign, n, rlen;
1683
1684 z = AriBuf;
1685 x = AriScratch;
1686 hilf = x + aribufSize;
1687
1688 n = bigretr(argStkPtr,x,&sign);
1689 if(n == aERROR) {
1690 error(isqrtsym,err_int,*argStkPtr);
1691 return(brkerr());
1692 }
1693 if(sign) {
1694 error(isqrtsym,err_p0num,*argStkPtr);
1695 return(brkerr());
1696 }
1697 n = bigsqrt(x,n,z,&rlen,hilf);
1698 return(mkint(0,z,n));
1699 }
1700 /*------------------------------------------------------------------*/
1701 /*
1702 ** returns 1 (resp. 0, -1) if number in *ptr1 is
1703 ** bigger than (resp. equal, smaller than) number in *ptr2
1704 */
cmpnums(ptr1,ptr2,type)1705 PUBLIC int cmpnums(ptr1,ptr2,type)
1706 truc *ptr1, *ptr2;
1707 int type; /* must be an integer or float type */
1708 {
1709 word2 *x, *y;
1710 int k, n1, n2, sign1, sign2, cmp;
1711
1712 if(*ptr1 == *ptr2)
1713 return(0);
1714 if(type == fFIXNUM) {
1715 sign1 = *SIGNPTR(ptr1);
1716 sign2 = *SIGNPTR(ptr2);
1717 if(sign1 != sign2)
1718 return(sign1 ? -1 : 1);
1719 cmp = (*WORD2PTR(ptr1) > *WORD2PTR(ptr2) ? 1 : -1);
1720 return(sign1 ? -cmp : cmp);
1721 }
1722 else if(type == fBIGNUM) {
1723 n1 = bigref(ptr1,&x,&sign1);
1724 n2 = bigref(ptr2,&y,&sign2);
1725 if(sign1 != sign2)
1726 return(sign1 ? -1 : 1);
1727 cmp = cmparr(x,n1,y,n2);
1728 return(sign1 ? -cmp : cmp);
1729 }
1730 else if(type == fGF2NINT) {
1731 n1 = bigref(ptr1,&x,&sign1);
1732 n2 = bigref(ptr2,&y,&sign2);
1733 return cmparr(x,n1,y,n2);
1734 }
1735 else if(type >= fFLTOBJ) {
1736 return(cmpfloats(ptr1,ptr2,fltprec(type)+1));
1737 }
1738 else
1739 return(aERROR);
1740 }
1741 /*-------------------------------------------------------------------*/
cmpfloats(ptr1,ptr2,prec)1742 PRIVATE int cmpfloats(ptr1,ptr2,prec)
1743 truc *ptr1, *ptr2;
1744 int prec;
1745 {
1746 numdata npt1, npt2;
1747 long expo1, expo2;
1748 int cmp, sign1, n1, n2;
1749
1750 npt1.digits = AriBuf;
1751 npt2.digits = AriScratch;
1752 getnumtrunc(prec,ptr1,&npt1);
1753 getnumtrunc(prec,ptr2,&npt2);
1754 sign1 = npt1.sign;
1755 if(sign1 != npt2.sign)
1756 return(sign1 ? -1 : 1);
1757 n1 = normfloat(prec,&npt1);
1758 n2 = normfloat(prec,&npt2);
1759 expo1 = npt1.expo;
1760 expo2 = npt2.expo;
1761 if(expo1 > expo2) {
1762 return(sign1 ? -1 : 1);
1763 }
1764 else if(expo1 < expo2) {
1765 return(sign1 ? 1 : -1);
1766 }
1767 else {
1768 cmp = cmparr(npt1.digits,n1,npt2.digits,n2);
1769 return(sign1 ? -cmp : cmp);
1770 }
1771 }
1772 /*-----------------------------------------------------------*/
1773 /*
1774 ** Num1 < Num2
1775 */
Farilt()1776 PRIVATE truc Farilt()
1777 {
1778 int cmp = Gcompare(ariltsym);
1779
1780 if(cmp == aERROR)
1781 return(brkerr());
1782 else
1783 return(cmp < 0 ? true : false);
1784 }
1785 /*-------------------------------------------------------------*/
1786 /*
1787 ** Num1 > Num2
1788 */
Farigt()1789 PRIVATE truc Farigt()
1790 {
1791 int cmp = Gcompare(arigtsym);
1792
1793 if(cmp == aERROR)
1794 return(brkerr());
1795 else
1796 return(cmp > 0 ? true : false);
1797 }
1798 /*-----------------------------------------------------------*/
1799 /*
1800 ** Num1 <= Num2
1801 */
Farile()1802 PRIVATE truc Farile()
1803 {
1804 int cmp = Gcompare(arilesym);
1805
1806 if(cmp == aERROR)
1807 return(brkerr());
1808 else
1809 return(cmp <= 0 ? true : false);
1810 }
1811 /*---------------------------------------------------------------*/
1812 /*
1813 ** Num1 >= Num2
1814 */
Farige()1815 PRIVATE truc Farige()
1816 {
1817 int cmp = Gcompare(arigesym);
1818
1819 if(cmp == aERROR)
1820 return(brkerr());
1821 else
1822 return(cmp >= 0 ? true : false);
1823 }
1824 /*---------------------------------------------------------------*/
Gcompare(symb)1825 PRIVATE int Gcompare(symb)
1826 truc symb;
1827 {
1828 truc obj;
1829 char *errmsg;
1830 char *str1, *str2;
1831 int type, type1;
1832
1833 type = *FLAGPTR(argStkPtr);
1834 if(type >= fFIXNUM) {
1835 if((type1 = *FLAGPTR(argStkPtr-1)) >= fFIXNUM) {
1836 type = (type >= type1 ? type : type1);
1837 return(cmpnums(argStkPtr-1,argStkPtr,type));
1838 }
1839 else {
1840 errmsg = err_mism;
1841 obj = argStkPtr[-1];
1842 }
1843 }
1844 else if(type == fSTRING) {
1845 if(*FLAGPTR(argStkPtr-1) == fSTRING) {
1846 str1 = STRINGPTR(argStkPtr-1);
1847 str2 = STRINGPTR(argStkPtr);
1848 return(strcmp(str1,str2));
1849 }
1850 else {
1851 errmsg = err_mism;
1852 obj = argStkPtr[-1];
1853 }
1854 }
1855 else if(type == fCHARACTER) {
1856 if(*FLAGPTR(argStkPtr-1) == fCHARACTER) {
1857 return(*WORD2PTR(argStkPtr-1) - *WORD2PTR(argStkPtr));
1858 }
1859 else {
1860 errmsg = err_mism;
1861 obj = argStkPtr[-1];
1862 }
1863 }
1864 else {
1865 errmsg = err_type;
1866 obj = argStkPtr[0];
1867 }
1868 return(error(symb,errmsg,obj));
1869 }
1870 /*--------------------------------------------------------*/
inirandstate(rr)1871 PRIVATE void inirandstate(rr)
1872 word2 *rr;
1873 {
1874 rr[1] = sysrand();
1875 nextrand(rr,3);
1876 rr[0] = sysrand();
1877 nextrand(rr,3);
1878 rr[3] = 1;
1879 }
1880 /*--------------------------------------------------------*/
nextrand(rr,n)1881 PRIVATE void nextrand(rr,n)
1882 word2 *rr;
1883 int n;
1884 {
1885 /* for compilers which don't understand 57777U */
1886 static unsigned inc = 57777, scal = 56857;
1887 /* 57777 = 1 mod 4, 56857 prime */
1888
1889 incarr(rr,n,inc);
1890 multarr(rr,n,scal,rr);
1891 rr[3] = 1;
1892 }
1893 /*------------------------------------------------------------------*/
1894 /*
1895 ** 2-byte random integer
1896 */
random2(u)1897 PUBLIC unsigned random2(u)
1898 unsigned u;
1899 {
1900 nextrand(RandSeed,2);
1901 return(RandSeed[1] % u);
1902 }
1903 /*------------------------------------------------------------------*/
1904 /*
1905 ** 4-byte random integer
1906 */
random4(u)1907 PUBLIC unsigned random4(u)
1908 unsigned u;
1909 {
1910 word4 v;
1911
1912 nextrand(RandSeed,3);
1913 v = big2long(RandSeed+1,2);
1914 return(v % u);
1915 }
1916 /*------------------------------------------------------------------*/
1917 /*
1918 ** random(bound)
1919 */
Frandom()1920 PRIVATE truc Frandom()
1921 {
1922 numdata acc, acc2;
1923 word2 *x;
1924 unsigned a, b;
1925 int i, n, m, prec, type;
1926
1927 type = chknum(randsym,argStkPtr);
1928 if(type == aERROR)
1929 return(brkerr());
1930 if(type == fFIXNUM) {
1931 nextrand(RandSeed,2);
1932 a = *WORD2PTR(argStkPtr);
1933 if(!a)
1934 return(zero);
1935 b = RandSeed[1] % a;
1936 return(mkfixnum(b));
1937 }
1938 else if(type >= fFLTOBJ) {
1939 curFltPrec = deffltprec();
1940 prec = curFltPrec + 1;
1941 for(x=AriBuf, i=0; i<curFltPrec; i+=2, x+=2) {
1942 nextrand(RandSeed,2);
1943 cpyarr(RandSeed,2,x);
1944 }
1945 acc.digits = AriBuf;
1946 acc.sign = 0;
1947 acc.len = curFltPrec;
1948 acc.expo = -(curFltPrec<<4);
1949 refnumtrunc(prec,argStkPtr,&acc2);
1950 n = multtrunc(prec,&acc,&acc2,AriScratch);
1951 return(mkfloat(curFltPrec,&acc));
1952 }
1953 else { /* bignum */
1954 n = *BIGLENPTR(argStkPtr);
1955 for(x=AriBuf, i=0; i<n; i+=2, x+=2) {
1956 nextrand(RandSeed,3);
1957 cpyarr(RandSeed+1,2,x);
1958 }
1959 m = n;
1960 while(m > 0 && AriBuf[m-1] == 0)
1961 m--;
1962 x = BIGNUMPTR(argStkPtr);
1963 m = modbig(AriBuf,m,x,n,AriScratch);
1964 return(mkint(0,AriBuf,m));
1965 }
1966 }
1967 /*------------------------------------------------------------------*/
Frandseed(argn)1968 PRIVATE truc Frandseed(argn)
1969 int argn;
1970 {
1971 word2 *x;
1972 int sign, n, m;
1973
1974 if(argn == 1) {
1975 n = bigref(argStkPtr,&x,&sign);
1976 if(n != aERROR) {
1977 m = (n > 3 ? 3 : n);
1978 cpyarr(x,m,RandSeed);
1979 setarr(RandSeed+m,3-m,0);
1980 }
1981 }
1982 return(mkint(0,RandSeed,4));
1983 }
1984 /*------------------------------------------------------------------*/
objfltprec(obj)1985 PRIVATE int objfltprec(obj)
1986 truc obj;
1987 {
1988 variant v;
1989 int flg, prec;
1990
1991 if(obj == sfloatsym || obj == dfloatsym ||
1992 obj == lfloatsym || obj == xfloatsym) {
1993 v.xx = SYMbind(obj);
1994 prec = v.pp.ww;
1995 }
1996 else {
1997 v.xx = obj;
1998 flg = v.pp.b0;
1999 if(flg >= fFLTOBJ)
2000 prec = fltprec(flg);
2001 else
2002 prec = deffltprec();
2003
2004 }
2005 return(prec);
2006 }
2007 /*------------------------------------------------------------------*/
Ffloat(argn)2008 PRIVATE truc Ffloat(argn)
2009 int argn; /* argn = 1 or 2 */
2010 {
2011 truc *argptr;
2012 numdata xx;
2013 int prec;
2014
2015 if(argn == 1)
2016 prec = deffltprec();
2017 else
2018 prec = precdesc(*argStkPtr);
2019
2020 argptr = argStkPtr - argn + 1;
2021 if(chknum(floatsym,argptr) == aERROR)
2022 return(*argptr);
2023
2024 xx.digits = AriBuf;
2025 getnumtrunc(prec+1,argptr,&xx);
2026
2027 return(mkfloat(prec,&xx));
2028 }
2029 /*------------------------------------------------------------------*/
precdesc(obj)2030 PRIVATE int precdesc(obj)
2031 truc obj;
2032 {
2033 variant v;
2034 int flg, prec, bits;
2035
2036 if(obj == sfloatsym || obj == dfloatsym ||
2037 obj == lfloatsym || obj == xfloatsym) {
2038 v.xx = SYMbind(obj);
2039 prec = v.pp.ww;
2040 }
2041 else {
2042 v.xx = obj;
2043 flg = v.pp.b0;
2044 if(flg == fFIXNUM) {
2045 bits = v.pp.ww;
2046 prec = (bits + 15)/16;
2047 }
2048 else
2049 prec = deffltprec();
2050 }
2051 return(prec);
2052 }
2053 /*------------------------------------------------------------------*/
Fsetfltprec(argn)2054 PRIVATE truc Fsetfltprec(argn)
2055 int argn;
2056 {
2057 int prec;
2058 truc obj;
2059 int prnflg = 1;
2060
2061 if(argn == 2) {
2062 obj = argStkPtr[-1];
2063 if(*argStkPtr == zero)
2064 prnflg = 0;
2065 }
2066 else
2067 obj = *argStkPtr;
2068 prec = precdesc(obj);
2069 prec = setfltprec(prec);
2070 if(prnflg)
2071 setprnprec(prec);
2072 return(mkfixnum(16*prec));
2073 }
2074 /*------------------------------------------------------------------*/
Fsetprnprec()2075 PRIVATE truc Fsetprnprec()
2076 {
2077 int prec;
2078
2079 prec = precdesc(*argStkPtr);
2080 prec = setprnprec(prec);
2081 return(mkfixnum(16*prec));
2082 }
2083 /*------------------------------------------------------------------*/
Fgetprnprec()2084 PRIVATE truc Fgetprnprec()
2085 {
2086 int prec;
2087
2088 prec = setprnprec(-1);
2089 return(mkfixnum(16*prec));
2090 }
2091 /*------------------------------------------------------------------*/
Fgetfltprec(argn)2092 PRIVATE truc Fgetfltprec(argn)
2093 int argn;
2094 {
2095 int prec;
2096
2097 if(argn == 0)
2098 prec = deffltprec();
2099 else
2100 prec = objfltprec(*argStkPtr);
2101
2102 return(mkfixnum(16*prec));
2103 }
2104 /*-------------------------------------------------------------*/
Fmaxfltprec()2105 PRIVATE truc Fmaxfltprec()
2106 {
2107 int prec;
2108
2109 prec = maxfltprec();
2110 return(mkfixnum(16*prec));
2111 }
2112 /*-------------------------------------------------------------*/
Fdecode()2113 PRIVATE truc Fdecode()
2114 {
2115 truc *ptr;
2116 truc vec, obj;
2117 numdata acc;
2118 long expo;
2119 int flg, len;
2120
2121 flg = *FLAGPTR(argStkPtr);
2122 if(flg < fFLTOBJ) {
2123 if(flg >= fFIXNUM) {
2124 *argStkPtr = Ffloat(1);
2125 flg = *FLAGPTR(argStkPtr);
2126 }
2127 else {
2128 error(decodsym,err_float,*argStkPtr);
2129 return(brkerr());
2130 }
2131 }
2132 len = fltprec(flg);
2133
2134 acc.digits = AriBuf;
2135 len = getnumtrunc(len,argStkPtr,&acc);
2136 expo = (len ? acc.expo : 0);
2137 obj = mkint(acc.sign,AriBuf,len);
2138 WORKpush(obj);
2139 obj = mkinum(expo);
2140 WORKpush(obj);
2141 vec = mkvect0(2);
2142 ptr = VECTOR(vec);
2143 ptr[1] = WORKretr();
2144 ptr[0] = WORKretr();
2145 return(vec);
2146 }
2147 /*-------------------------------------------------------------*/
Fbitnot()2148 PRIVATE truc Fbitnot()
2149 {
2150 int n, sign;
2151
2152 if(chkints(bnotsym,argStkPtr,1) == aERROR)
2153 return(brkerr());
2154
2155 n = bigretr(argStkPtr,AriBuf,&sign);
2156 if(sign) {
2157 n = decarr(AriBuf,n,1);
2158 sign = 0;
2159 }
2160 else {
2161 n = incarr(AriBuf,n,1);
2162 sign = 1;
2163 }
2164 return(mkint(sign,AriBuf,n));
2165 }
2166 /*-------------------------------------------------------------*/
Fbitset()2167 PRIVATE truc Fbitset()
2168 {
2169 return(Gbitset(bsetsym));
2170 }
2171 /*-------------------------------------------------------------*/
Fbitclear()2172 PRIVATE truc Fbitclear()
2173 {
2174 return(Gbitset(bclrsym));
2175 }
2176 /*-------------------------------------------------------------*/
Gbitset(symb)2177 PRIVATE truc Gbitset(symb)
2178 truc symb;
2179 {
2180 word2 *x, *y;
2181 long index;
2182 int b, n, n1, m, sign, sign1;
2183 word2 u, mask = 1;
2184
2185 if(chkints(symb,argStkPtr-1,2) == aERROR)
2186 return(brkerr());
2187 x = AriBuf;
2188 n = twocretr(argStkPtr-1,x);
2189 u = x[n];
2190 sign = (u == 0xFFFF ? MINUSBYTE : 0);
2191
2192 n1 = bigref(argStkPtr,&y,&sign1);
2193 if(sign1) {
2194 error(symb,err_p0num,*argStkPtr);
2195 return(brkerr());
2196 }
2197 if(n1 > 2)
2198 index = 0xFFFF0;
2199 else
2200 index = big2long(y,n1);
2201 b = index & 0xF;
2202 index >>= 4;
2203 m = index;
2204 if(index >= n) {
2205 if((sign && symb == bsetsym) ||
2206 (!sign && symb == bclrsym)) {
2207 return(argStkPtr[-1]); /* no action taken */
2208 }
2209 else if(index >= aribufSize) {
2210 error(symb,err_ovfl,voidsym);
2211 return(argStkPtr[-1]);
2212 }
2213 else {
2214 setarr(x+n+1,m-n,u);
2215 n = m + 1;
2216 }
2217 }
2218 mask <<= b;
2219 if(symb == bsetsym)
2220 x[m] |= mask;
2221 else
2222 x[m] &= ~mask;
2223 while((n > 0) && (x[n-1] == u))
2224 n--;
2225 if(sign) {
2226 notarr(x,n);
2227 n = incarr(x,n,1);
2228 }
2229 return(mkint(sign,x,n));
2230 }
2231 /*-------------------------------------------------------------*/
Fbittest()2232 PRIVATE truc Fbittest()
2233 {
2234 word2 *x;
2235 long index;
2236 int k, b, n, sign;
2237 word2 mask = 1;
2238
2239 if(chkintt(btestsym,argStkPtr-1) == aERROR ||
2240 chkint(btestsym,argStkPtr) == aERROR)
2241 return(brkerr());
2242
2243 n = bigref(argStkPtr-1,&x,&sign);
2244 if(sign) {
2245 x = AriBuf;
2246 n = twocretr(argStkPtr-1,x);
2247 }
2248 index = intretr(argStkPtr);
2249 if(index == LONGERROR)
2250 return(sign ? constone : zero);
2251 else if(index < 0) {
2252 error(btestsym,err_p0num,*argStkPtr);
2253 return(brkerr());
2254 }
2255 k = index >> 4;
2256 if(index > 0x7FFF0 || k >= n)
2257 return(sign ? constone : zero);
2258 b = index & 0xF;
2259 mask <<= b;
2260 return(x[k] & mask ? constone : zero);
2261 }
2262 /*-------------------------------------------------------------*/
Fbitshift()2263 PRIVATE truc Fbitshift()
2264 {
2265 long sh, nn;
2266 int n, sign;
2267
2268 if(chkints(bshiftsym,argStkPtr-1,2) == aERROR)
2269 return(brkerr());
2270
2271 nn = n = bigretr(argStkPtr-1,AriBuf,&sign);
2272 sh = intretr(argStkPtr);
2273 if(sh == LONGERROR || sh >= maxfltex - (nn<<4)) {
2274 error(bshiftsym,err_ovfl,voidsym);
2275 return(brkerr());
2276 }
2277 if(sign && sh < 0) {
2278 n = decarr(AriBuf,n,1);
2279 n = lshiftarr(AriBuf,n,sh);
2280 n = incarr(AriBuf,n,1);
2281 }
2282 else
2283 n = lshiftarr(AriBuf,n,sh);
2284 return(mkint(sign,AriBuf,n));
2285 }
2286 /*-------------------------------------------------------------*/
Fbitlength()2287 PRIVATE truc Fbitlength()
2288 {
2289 word2 *x;
2290 long len;
2291 int n, sign;
2292
2293 if(chkintt(blensym,argStkPtr) == aERROR)
2294 return(brkerr());
2295
2296 n = bigref(argStkPtr,&x,&sign);
2297 len = bit_length(x,n);
2298 return(mkinum(len));
2299 }
2300 /*-------------------------------------------------------------*/
Fbitcount()2301 PRIVATE truc Fbitcount()
2302 {
2303 word2 *x;
2304 long count;
2305 int k, n, sign;
2306
2307 if(chkintt(blensym,argStkPtr) == aERROR)
2308 return(brkerr());
2309
2310 n = bigref(argStkPtr,&x,&sign);
2311 count = 0;
2312 for(k=0; k<n; k++)
2313 count += bitcount(x[k]);
2314 return(mkinum(count));
2315 }
2316 /*-------------------------------------------------------------*/
Fcardinal()2317 PRIVATE truc Fcardinal()
2318 {
2319 return(Gcardinal(cardsym));
2320 }
2321 /*-------------------------------------------------------------*/
Finteger()2322 PRIVATE truc Finteger()
2323 {
2324 return(Gcardinal(int_sym));
2325 }
2326 /*-------------------------------------------------------------*/
Gcardinal(symb)2327 PRIVATE truc Gcardinal(symb)
2328 truc symb;
2329 {
2330 word2 *x;
2331 byte *bpt;
2332 unsigned u;
2333 unsigned len;
2334 int i, n, flg;
2335 int sign = 0;
2336
2337 flg = *FLAGPTR(argStkPtr);
2338 if(flg == fGF2NINT) {
2339 n = bigretr(argStkPtr,AriBuf,&sign);
2340 return mkint(0,AriBuf,n);
2341 }
2342 if(flg != fBYTESTRING && flg != fSTRING) {
2343 error(symb,err_bystr,*argStkPtr);
2344 return(brkerr());
2345 }
2346 len = *STRLENPTR(argStkPtr);
2347 if(len >= aribufSize*2 - 2) {
2348 error(symb,err_2long,mkfixnum(len));
2349 return(brkerr());
2350 }
2351 bpt = (byte *)STRINGPTR(argStkPtr);
2352 if(symb == int_sym && (bpt[len-1] & 0x80))
2353 sign = MINUSBYTE;
2354 n = len / 2;
2355 x = AriBuf;
2356 for(i=0; i<n; i++) {
2357 u = bpt[1];
2358 u <<= 8;
2359 u += bpt[0];
2360 *x++ = u;
2361 bpt += 2;
2362 }
2363 if(len & 1) {
2364 *x = *bpt;
2365 if(sign)
2366 *x |= 0xFF00;
2367 n++;
2368 }
2369 if(sign) {
2370 notarr(AriBuf,n);
2371 n = incarr(AriBuf,n,1);
2372 }
2373 while(n > 0 && AriBuf[n-1] == 0)
2374 n--;
2375 return(mkint(sign,AriBuf,n));
2376 }
2377 /*-------------------------------------------------------------*/
Fbitand()2378 PRIVATE truc Fbitand()
2379 {
2380 return(Gboole(bandsym,and2arr));
2381 }
2382 /*-------------------------------------------------------------*/
Fbitor()2383 PRIVATE truc Fbitor()
2384 {
2385 return(Gboole(borsym,or2arr));
2386 }
2387 /*-------------------------------------------------------------*/
Fbitxor()2388 PRIVATE truc Fbitxor()
2389 {
2390 return(Gboole(bxorsym,xor2arr));
2391 }
2392 /*------------------------------------------------------------------*/
Gboole(symb,boolfun)2393 PRIVATE truc Gboole(symb,boolfun)
2394 truc symb;
2395 ifunaa boolfun;
2396 {
2397 int sign;
2398 int n, m;
2399
2400 if(chkints(symb,argStkPtr-1,2) == aERROR)
2401 return(brkerr());
2402 n = twocretr(argStkPtr-1,AriBuf);
2403 m = twocretr(argStkPtr,AriScratch);
2404 n = boolfun(AriBuf,n,AriScratch,m);
2405
2406 sign = (AriBuf[n] == 0xFFFF ? MINUSBYTE : 0);
2407 if(sign) {
2408 notarr(AriBuf,n);
2409 n = incarr(AriBuf,n,1);
2410 }
2411 return(mkint(sign,AriBuf,n));
2412 }
2413 /*--------------------------------------------------------*/
chkplusargs(sym,argptr)2414 PRIVATE int chkplusargs(sym,argptr)
2415 truc sym;
2416 truc *argptr;
2417 {
2418 int flg, flg1;
2419
2420 flg = *FLAGPTR(argptr);
2421 flg1 = *FLAGPTR(argptr+1);
2422 if(flg >= fFIXNUM) {
2423 if(flg1 >= fFIXNUM)
2424 return (flg1 >= flg ? flg1 : flg);
2425 else
2426 return error(sym,err_num,argptr[1]);
2427 }
2428 else if((flg == flg1) && (flg == fVECTOR || flg == fGF2NINT))
2429 return flg;
2430 else
2431 return error(sym,err_num,*argptr);
2432 }
2433 /*--------------------------------------------------------*/
chktimesargs(argptr)2434 PRIVATE int chktimesargs(argptr)
2435 truc *argptr;
2436 {
2437 int flg, flg1;
2438 truc ele;
2439
2440 flg = *FLAGPTR(argptr);
2441 flg1 = *FLAGPTR(argptr+1);
2442
2443 if(flg < fFIXNUM) {
2444 if(flg == fVECTOR) { /* then second argument must be a scalar */
2445 if(flg1 >= fFIXNUM || flg1 == fGF2NINT) { /* swap args */
2446 ele = *argptr;
2447 *argptr = argptr[1];
2448 *(argptr+1) = ele;
2449 return (flg1 | (flg << 8));
2450 }
2451 else
2452 return error(timessym,err_num,argptr[1]);
2453 }
2454 else if(flg == fGF2NINT) {
2455 if(flg1 == flg)
2456 return flg;
2457 else if(flg1 == fVECTOR)
2458 return (flg | (flg1 << 8));
2459 }
2460 return error(timessym,err_num,*argptr);
2461 }
2462 /* here flg >= fFIXNUM */
2463 if(flg1 >= fFIXNUM) {
2464 return (flg1 >= flg ? flg1 : flg);
2465 }
2466 else if(flg1 == fVECTOR)
2467 return (flg | (flg1 << 8));
2468 else
2469 return error(timessym,err_num,argptr[1]);
2470 }
2471 /*--------------------------------------------------------*/
chkmodargs(sym,argptr)2472 PRIVATE int chkmodargs(sym,argptr)
2473 truc sym;
2474 truc *argptr;
2475 {
2476 int flg, flg0;
2477
2478 flg0 = *FLAGPTR(argptr+1);
2479 if(flg0 < fFIXNUM || flg0 > fBIGNUM)
2480 return error(sym,err_num,argptr[1]);
2481 flg = *FLAGPTR(argptr);
2482 if((flg < fFIXNUM && flg != fVECTOR) || (flg > fBIGNUM))
2483 return error(sym,err_num,argptr[0]);
2484 if(flg == fVECTOR)
2485 return((fVECTOR<<8) | flg0);
2486 else
2487 return(flg0);
2488 }
2489 /*--------------------------------------------------------*/
chkdivfargs(sym,argptr)2490 PRIVATE int chkdivfargs(sym,argptr)
2491 truc sym;
2492 truc *argptr;
2493 {
2494 int flg, flg0;
2495
2496 flg0 = *FLAGPTR(argptr+1);
2497 flg = *FLAGPTR(argptr);
2498 if(flg0 < fFIXNUM) {
2499 if(flg0 == fGF2NINT && flg == flg)
2500 return flg;
2501 else
2502 return error(sym,err_num,argptr[1]);
2503 }
2504 if(flg < fFIXNUM && flg != fVECTOR)
2505 return error(sym,err_num,argptr[0]);
2506
2507 if(flg == fVECTOR)
2508 return((fVECTOR<<8) | flg0);
2509 else
2510 return(flg0);
2511 }
2512 /*********************************************************************/
2513