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,&quot);
1350     m = refnumtrunc(prec,ptr+1,&temp);
1351     n = divtrunc(prec,&quot,&temp,AriScratch);
1352     if(n < 0) {
1353         error(divfsym,err_ovfl,voidsym);
1354         return(brkerr());
1355     }
1356     return(mkfloat(curFltPrec,&quot));
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