123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692 |
- /* arith03.c Copyright (C) 1990-96 Codemist Ltd */
- /*
- * Arithmetic functions.
- * quotient.
- *
- */
- /* Signature: 7f5d27b7 31-May-1997 */
- #include <stdarg.h>
- #include <string.h>
- #include <ctype.h>
- #include <math.h>
- #include "machine.h"
- #include "tags.h"
- #include "cslerror.h"
- #include "externs.h"
- #include "arith.h"
- #ifdef TIMEOUT
- #include "timeout.h"
- #endif
- /*
- * Division
- */
- #ifndef IDIVIDE
- #ifdef MULDIV64
- unsigned32 Idivide(unsigned32 *qp, unsigned32 a, unsigned32 b, unsigned32 c)
- /*
- * *qp = (a,b) / c, return the remainder
- *
- * The double-length value (a,b) is given as a 62-bit positive number. Its
- * low order 31 bits are in b, and the 0x80000000 bit of b is zero. The
- * high order part (in a) is also represented with the 0x80000000 bit zero.
- * The divisor c is a positive value that is at most 0x7fffffff, and the
- * procedure will only be called when the correct quotient is in the
- * (inclusive) range 0 to 0x7fffffff. The above constraints can be thought
- * of in two ways: (a) Idivide used a 31-bit word and is working on
- * unsigned values. The 31-bit quantities happen to be being passed around
- * in 32 bit registers, but the top bit of those registers will never be used
- * and will contain zero, or (b) the range of values used is such that
- * a 64 by 32-bit division can be performed, and by constraining the range
- * of values this 64 by 32-bit division only has to handle positive numbers,
- * but depending on the hardware of your computer you are entitled to use
- * either signed or unsigned arithmetic.
- */
- {
- unsigned64 p = ((unsigned64)a << 31) | (unsigned64)b;
- *qp = (unsigned32)(p / (unsigned64)c);
- return (unsigned32)(p % (unsigned64)c);
- }
- #else
- unsigned32 Idivide(unsigned32 *qp, unsigned32 a, unsigned32 b, unsigned32 c)
- /*
- * *qp = (a,b) / c, return the remainder
- *
- * The double-length value (a,b) is given as a 62-bit positive number. Its
- * low order 31 bits are in b, and the 0x80000000 bit of b is zero. The
- * high order part (in a) is also represented with the 0x80000000 bit zero.
- * The divisor c is a positive value that is at most 0x7fffffff, and the
- * procedure will only be called when the correct quotient is in the
- * (inclusive) range 0 to 0x7fffffff. The above constraints can be thought
- * of in two ways: (a) Idivide used a 31-bit word and is working on
- * unsigned values. The 31-bit quantities happen to be being passed around
- * in 32 bit registers, but the top bit of those registers will never be used
- * and will contain zero, or (b) the range of values used is such that
- * a 64 by 32-bit division can be performed, and by constraining the range
- * of values this 64 by 32-bit division only has to handle positive numbers,
- * but depending on the hardware of your computer you are entitled to use
- * either signed or unsigned arithmetic.
- */
- {
- int i;
- unsigned32 q = 0;
- /*
- * I have a loop that needs to be executed 32 times, with just
- * slightly different behaviour the last time around. Since it is
- * fairly short, I unwind it three-fold into a loop executed 10 times
- * plus a postlude. I also do a test at the start that decides if the
- * quotient will be small (less than about 16 bits) an in that case
- * loop rather less often - my reasoning is that a plausible distribution
- * of quotients is exponential so the short route will be taken about
- * half of the time, and will save almost half of the work done here at the
- * cost of a not-too-expensive extra test to begin with.
- */
- if (a < (c >> 15))
- { a = (a << 15) | (b >> 16);
- b = (b << 15) & 0x7fffffff;
- i = 5;
- }
- else i = 0;
- do
- { q = q << 3; /* accumulate quotient 3 bits at a time */
- if (a >= c) /* bit 1 of 3... */
- { a = a - c;
- q = q + 4;
- }
- a = a << 1; /* shift (a,b) doubleword left 1 bit */
- b = b << 1;
- if (((int32)b) < 0) a = a + 1;
- if (a >= c) /* bit 2 of 3... */
- { a = a - c;
- q = q + 2;
- }
- a = a << 1;
- b = b << 1;
- if (((int32)b) < 0) a = a + 1;
- if (a >= c) /* bit 3 of 3... */
- { a = a - c;
- q = q + 1;
- }
- a = a << 1;
- b = b << 1;
- if (((int32)b) < 0) a = a + 1;
- i++;
- } while (i < 10);
- q = q << 2; /* 2 more bits to be included */
- if (a >= c)
- { a = a - c;
- q = q + 2;
- }
- a = a << 1;
- b = b << 1;
- if (((int32)b) < 0) a = a + 1;
- if (a >= c) /* the final bit of the quotient */
- { a = a - c; /* leave remainder in a, b should be zero */
- q = q + 1;
- }
- *qp = q;
- return a;
- }
- #endif
- #endif /* IDIVIDE */
- #ifdef COMMON
- static Lisp_Object quotis(Lisp_Object a, Lisp_Object b)
- {
- Float_union bb;
- bb.i = b - TAG_SFLOAT;
- if (bb.f == 0.0) return aerror2("bad arg for quotient", a, b);
- bb.f = (float) ((float)int_of_fixnum(a) / bb.f);
- return (bb.i & ~(int32)0xf) + TAG_SFLOAT;
- }
- #endif
- static Lisp_Object quotib(Lisp_Object a, Lisp_Object b)
- {
- /*
- * a fixnum divided by a bignum always gives 0, regardless of signs etc.,
- * save in the one case of (-#x8000000)/(#x8000000) which must yield -1
- */
- if (int_of_fixnum(a) == fix_mask && bignum_length(b) == 8 &&
- bignum_digits(b)[0] == 0x08000000) return fixnum_of_int(-1);
- else return fixnum_of_int(0);
- }
- #ifdef COMMON
- static Lisp_Object CLquotib(Lisp_Object a, Lisp_Object b)
- {
- Lisp_Object g, nil = C_nil;
- bool w;
- push2(a, b);
- w = minusp(b);
- errexitn(2);
- g = gcd(stack[0], stack[-1]);
- errexitn(2);
- if (w)
- { g = negate(g);
- errexitn(2);
- }
- a = stack[-1];
- push(g);
- a = quot2(a, g);
- errexitn(3);
- pop2(g, b);
- stack[0] = a;
- b = quot2(b, g);
- pop(a);
- errexit();
- return make_ratio(a, b);
- }
- static Lisp_Object CLquotbi(Lisp_Object a, Lisp_Object b)
- {
- return CLquotib(a, b);
- }
- static Lisp_Object CLquotbb(Lisp_Object a, Lisp_Object b)
- {
- return CLquotib(a, b);
- }
- static Lisp_Object quotir(Lisp_Object a, Lisp_Object b)
- {
- Lisp_Object w, nil;
- if (a == fixnum_of_int(0)) return a;
- push3(b, a, C_nil);
- #define g stack[0]
- #define a stack[-1]
- #define b stack[-2]
- g = gcd(a, numerator(b));
- /*
- * the gcd() function guarantees to hand back a positive result.
- */
- nil = C_nil;
- if (exception_pending()) goto fail;
- a = quot2(a, g);
- nil = C_nil;
- if (exception_pending()) goto fail;
- w = minusp(numerator(b));
- nil = C_nil;
- if (exception_pending()) goto fail;
- if (w)
- { g = negate(g);
- nil = C_nil;
- if (exception_pending()) goto fail;
- }
- g = quot2(numerator(b), g); /* denominator of result will be +ve */
- nil = C_nil;
- if (exception_pending()) goto fail;
- a = times2(a, denominator(b));
- nil = C_nil;
- if (exception_pending()) goto fail;
- w = make_ratio(a, g);
- popv(3);
- return w;
- fail:
- popv(3);
- return nil;
- #undef a
- #undef b
- #undef g
- }
- static Lisp_Object quotic(Lisp_Object a, Lisp_Object b)
- /*
- * Used for all cases of xxx/<p+iq>. This is coded in a fairly naive
- * way, multiplying both numerator and denominator of what will end up
- * as the result by the complex conjugate of b. If floating point
- * arithmetic is used this can lead to grossly premature overflow. For
- * /* the moment I will ignore that miserable fact
- */
- {
- Lisp_Object u, v, nil;
- push2(a, b);
- #define b stack[0]
- #define a stack[-1]
- /*
- * a / (p + iq) is computed as follows:
- * (a * (p - iq)) / (p^2 + q^2)
- */
- u = negate(imag_part(b));
- nil = C_nil;
- if (exception_pending()) goto fail;
- u = make_complex(real_part(b), u);
- nil = C_nil;
- if (exception_pending()) goto fail;
- a = times2(a, u);
- u = real_part(b);
- u = times2(u, u);
- nil = C_nil;
- if (exception_pending()) goto fail;
- v = imag_part(b);
- b = u;
- u = times2(v, v);
- nil = C_nil;
- if (exception_pending()) goto fail;
- u = plus2(u, b);
- nil = C_nil;
- if (exception_pending()) goto fail;
- v = a;
- popv(2);
- return quot2(v, u);
- #undef a
- #undef b
- fail:
- popv(2);
- return nil;
- }
- #endif
- static Lisp_Object quotif(Lisp_Object a, Lisp_Object b)
- {
- double d = float_of_number(b);
- if (d == 0.0) return aerror2("bad arg for quotient", a, b);
- d = (double)int_of_fixnum(a) / d;
- return make_boxfloat(d, type_of_header(flthdr(b)));
- }
- #ifdef COMMON
- static Lisp_Object quotsi(Lisp_Object a, Lisp_Object b)
- {
- Float_union aa;
- if (b == fixnum_of_int(0)) return aerror2("bad arg for quotient", a, b);
- aa.i = a - TAG_SFLOAT;
- aa.f = (float) (aa.f / (float)int_of_fixnum(b));
- return (aa.i & ~(int32)0xf) + TAG_SFLOAT;
- }
- static Lisp_Object quotsb(Lisp_Object a, Lisp_Object b)
- {
- double d = float_of_number(b);
- if (d == 0.0) return aerror2("bad arg for quotient", a, b);
- d = float_of_number(a) / d;
- return make_sfloat(d);
- }
- #define quotsr(a, b) quotsb(a, b)
- #define quotsc(a, b) quotic(a, b)
- #endif
- static Lisp_Object quotsf(Lisp_Object a, Lisp_Object b)
- {
- double d = float_of_number(b);
- if (d == 0.0) return aerror2("bad arg for quotient", a, b);
- d = float_of_number(a) / d;
- return make_boxfloat(d, type_of_header(flthdr(b)));
- }
- Lisp_Object quotbn(Lisp_Object a, int32 n)
- /*
- * Divide a bignum by an integer, where the integer is (by now)
- * a natural C int32 but limited to 31 not 32 bits active. I.e.
- * this can be used when dividing by a fixnum or by dividing by
- * a one-word bignum. I will not call this with n=-1, which would
- * be the one case where it could cause the quotient to be bigger
- * than the dividend. Leaves nwork set to the remainder.
- */
- {
- Lisp_Object nil;
- int sign;
- int32 lena = (bignum_length(a)>>2)-2, i, lenc, lenx;
- unsigned32 carry;
- if (lena == 0) /* one-word bignum as numerator */
- { int32 p = (int32)bignum_digits(a)[0], w;
- nil = C_nil; /* So I can access nwork */
- nwork = p % n;
- /*
- * C does not define what happens on non-exact division involving
- * negative quantities, so I adjust things here so that the remainder
- * has the sign that I want and the division I do is exact.
- */
- if (p < 0)
- { if (nwork > 0) nwork -= n;
- }
- else if (nwork < 0) nwork += n;
- /* Remainder should be correct (with correct sign) by now, regardless */
- p = (p - nwork) / n;
- w = p & fix_mask;
- if (w == 0 || w == fix_mask) return fixnum_of_int(p);
- else return make_one_word_bignum(p);
- }
- else if (lena == 1)
- /*
- * I treat division of a 2-word bignum by a fixnum or 1-word bignum as
- * a special case since it can lead to a fixnum result - if the divisor is
- * just one word long and the dividend is 3 or more words I would
- * certainly have a bignum result. Thus by separating off the code here
- * I (a) remove the need for test relating to big- to fixnum conversion
- * later on and (b) avoid allocating heap in this tolerably common case
- * when sometimes I will not need to.
- */
- { unsigned32 a0 = bignum_digits(a)[0];
- int32 a1 = (int32)bignum_digits(a)[1];
- if (a1 < 0)
- { sign = 3;
- if (a0 == 0) a1 = -a1;
- else
- { a0 = clear_top_bit(-(int32)a0);
- a1 = ~a1;
- }
- }
- else sign = 0;
- if (n < 0) sign ^= 1, n = -n;
- if (a1 < n) /* result can be calculated in one word - good! */
- { unsigned32 q, r;
- Ddivide(r, q, (unsigned32)a1, a0, n);
- if ((sign & 2) != 0) r = -(int32)r;
- nil = C_nil;
- nwork = r;
- a0 = q;
- if (a0 == 0) return fixnum_of_int(0);
- if ((sign & 1) == 0)
- { if ((a0 & fix_mask) == 0) return fixnum_of_int(a0);
- else if ((a0 & 0x40000000) == 0)
- return make_one_word_bignum(a0);
- else return make_two_word_bignum(0, a0);
- }
- a0 = -(int32)a0;
- if ((a0 & fix_mask) == fix_mask) return fixnum_of_int(a0);
- else if ((a0 & 0xc0000000U) == 0xc0000000U)
- return make_one_word_bignum(a0);
- else return make_two_word_bignum(-1, clear_top_bit(a0));
- }
- /*
- * here the quotient will involve two digits.
- */
- else
- { unsigned32 q1, q0, r;
- Ddivide(r, q1, 0, (unsigned32)a1, n);
- Ddivide(r, q0, r, a0, n);
- if ((sign & 2) != 0) r = -(int32)r;
- nil = C_nil;
- nwork = r;
- if ((sign & 1) != 0)
- { if (q0 == 0) q1 = -(int32)q1;
- else
- { q1 = ~q1;
- q0 = clear_top_bit(-(int32)q0);
- }
- }
- if (q1 == 0 && (q0 & 0x40000000) == 0)
- return make_one_word_bignum(q0);
- return make_two_word_bignum(q1, q0);
- }
- }
- /*
- * That has dealt with the special cases - now the input is a bignum
- * with at least 3 digits and the quotient will certainly be a bignum.
- * Start by allocating a workspace copy of the dividend. negateb will
- * leave a a bignum, although it may change its length.
- */
- if ((int32)bignum_digits(a)[lena] < 0) a = negateb(a), sign = 3;
- else a = copyb(a), sign = 0;
- errexit();
- if (n < 0)
- { sign ^= 1;
- n = -n;
- }
- /*
- * Now both a and n have been made positive, and their original signs
- * have been recorded so that I can tidy up at the end. The 1 bit in sign
- * tells me what sign the quotient will have, the 2 bit gives the sign
- * for the remainder.
- */
- lena = (bignum_length(a)>>2)-1;
- carry = 0;
- for (i=lena-1; i>=0; i--)
- Ddivide(carry, bignum_digits(a)[i], carry, bignum_digits(a)[i], n);
- if ((sign & 2) != 0) carry = -(int32)carry;
- nil = C_nil;
- nwork = carry; /* leave remainder available to caller */
- lena--;
- while (bignum_digits(a)[lena] == 0) lena--;
- if ((bignum_digits(a)[lena] & 0x40000000) != 0) lena++;
- if ((sign & 1) != 0) /* quotient will be negative */
- { carry = 0xffffffffU;
- for (i=0; i<lena; i++)
- { carry = clear_top_bit(~bignum_digits(a)[i]) + top_bit(carry);
- bignum_digits(a)[i] = clear_top_bit(carry);
- }
- carry = ~bignum_digits(a)[i] + top_bit(carry);
- if (carry == -1 && (bignum_digits(a)[i-1] & 0x40000000) != 0)
- { bignum_digits(a)[lena] = 0;
- lena--;
- bignum_digits(a)[i-1] |= ~0x7fffffff;
- }
- else bignum_digits(a)[i] = carry;
- }
- /* I need to free up some space here, I guess */
- lenx = (bignum_length(a)>>2)-1;
- lenx |= 1; /* round up to allow for allocation in doublewords */
- lenc = (lena+1) | 1;
- if (lenc != lenx) /* space to discard? */
- bignum_digits(a)[lenc] = make_bighdr(lenx-lenc);
- numhdr(a) = make_bighdr(lena+2);
- return a;
- }
- Lisp_Object quotbn1(Lisp_Object a, int32 n)
- /*
- * Divide a bignum by an integer, where the integer is (by now)
- * a natural C int32 but limited to 31 not 32 bits active. I.e.
- * this can be used when dividing by a fixnum or by dividing by
- * a one-word bignum. I will not call this with n=-1, which would
- * be the one case where it could cause the quotient to be bigger
- * than the dividend. Leaves nwork set to the remainder.
- * This version is JUST the same as quotbn() except that it does not
- * hand back a useful quotient - it is intended for use when only
- * the remainder is wanted. For consistency I leave that where quotbn() did.
- */
- {
- Lisp_Object nil;
- int sign;
- int32 lena = (bignum_length(a)>>2)-2, i;
- unsigned32 carry;
- if (lena == 0) /* one-word bignum as numerator */
- { int32 p = (int32)bignum_digits(a)[0];
- nil = C_nil; /* So I can access nwork */
- nwork = p % n;
- /*
- * C does not define what happens on non-exact division involving
- * negative quantities, so I adjust things here so that the remainder
- * has the sign that I want and the division I do is exact.
- */
- if (p < 0)
- { if (nwork > 0) nwork -= n;
- }
- else if (nwork < 0) nwork += n;
- /* Remainder should be correct (with correct sign) by now, regardless */
- else return nil;
- }
- else if (lena == 1)
- /*
- * I treat division of a 2-word bignum by a fixnum or 1-word bignum as
- * a special case since it can lead to a fixnum result - if the divisor is
- * just one word long and the dividend is 3 or more words I would
- * certainly have a bignum result. Thus by separating off the code here
- * I (a) remove the need for test relating to big- to fixnum conversion
- * later on and (b) avoid allocating heap in this tolerably common case
- * when sometimes I will not need to.
- */
- { unsigned32 a0 = bignum_digits(a)[0];
- int32 a1 = (int32)bignum_digits(a)[1];
- if (a1 < 0)
- { sign = 3;
- if (a0 == 0) a1 = -a1;
- else
- { a0 = clear_top_bit(-(int32)a0);
- a1 = ~a1;
- }
- }
- else sign = 0;
- if (n < 0) sign ^= 1, n = -n;
- if (a1 < n) /* result can be calculated in one word - good! */
- { unsigned32 q, r;
- Ddivide(r, q, (unsigned32)a1, a0, n);
- if ((sign & 2) != 0) r = -(int32)r;
- nil = C_nil;
- nwork = r;
- return nil;
- }
- /*
- * here the quotient will involve two digits.
- */
- else
- { unsigned32 q1, q0, r;
- Ddivide(r, q1, 0, (unsigned32)a1, n);
- Ddivide(r, q0, r, a0, n);
- if ((sign & 2) != 0) r = -(int32)r;
- nil = C_nil;
- nwork = r;
- return nil;
- }
- }
- /*
- * That has dealt with the special cases - now the input is a bignum
- * with at least 3 digits and the quotient will certainly be a bignum.
- * Start by allocating a workspace copy of the dividend. negateb will
- * leave a a bignum, although it may change its length.
- */
- /*
- * Also note that I could fold the negateb() in with the division, and
- * thereby save allocating a big hunk of extra memory.
- */
- if ((int32)bignum_digits(a)[lena] < 0) a = negateb(a), sign = 3;
- else a = copyb(a), sign = 0;
- errexit();
- if (n < 0)
- { sign ^= 1;
- n = -n;
- }
- /*
- * Now both a and n have been made positive, and their original signs
- * have been recorded so that I can tidy up at the end. The 1 bit in sign
- * tells me what sign the quotient will have, the 2 bit gives the sign
- * for the remainder.
- */
- lena = (bignum_length(a)>>2)-1;
- carry = 0;
- for (i=lena-1; i>=0; i--)
- Ddivide(carry, bignum_digits(a)[i], carry, bignum_digits(a)[i], n);
- if ((sign & 2) != 0) carry = -(int32)carry;
- nil = C_nil;
- nwork = carry; /* leave remainder available to caller */
- return nil;
- }
- static Lisp_Object quotbi(Lisp_Object a, Lisp_Object b)
- {
- /*
- * dividing by 1 or -1 seems worth optimising.
- */
- if (b == fixnum_of_int(1)) return a;
- else if (b == fixnum_of_int(-1)) return negateb(a);
- else if (b == fixnum_of_int(0))
- return aerror2("bad arg for quotient", a, b);
- else return quotbn(a, int_of_fixnum(b));
- }
- #ifdef COMMON
- #define quotbs(a, b) quotsb(a, b)
- #endif
- #ifdef __powerc
- /* If you have trouble compiling this just comment it out, please */
- #pragma options(!global_optimizer)
- #endif
- Lisp_Object quotbb(Lisp_Object a, Lisp_Object b)
- /*
- * Divide one bignum by another. This can compute the
- * remainder at the same time as the quotient, and leaves same around
- * in mv_2. If it is not needed then setting up the remainder is
- * a cost - but usually a modest one (in context), so I think that
- * the simplification is worth-while.
- */
- {
- Lisp_Object nil, olda;
- int sign;
- int32 lena, lenb, lenc, lenx, i;
- unsigned32 carry, scale, v1;
- Lisp_Object q, r;
- /*
- * If I am dividing by a one-word bignum I can go a bit faster...
- */
- lenb = (bignum_length(b)>>2)-2;
- if (lenb == 0)
- { q = quotbn(a, bignum_digits(b)[0]);
- errexit();
- /*
- * Now lots of frivolity packing up the remainder...
- */
- r = (Lisp_Object)(nwork & fix_mask);
- if ((int32)r == 0 || (int32)r == fix_mask)
- mv_2 = fixnum_of_int(nwork);
- else
- { push(q);
- a = make_one_word_bignum(nwork);
- pop(q);
- errexit();
- mv_2 = a;
- }
- return q;
- }
- /*
- * Convert to sign and magnitude representation
- */
- push2(a, b);
- lena = (bignum_length(a)>>2)-2;
- if ((int32)bignum_digits(a)[lena] < 0)
- { a = negateb(a);
- sign = 3;
- }
- else
- { a = copyb(a);
- sign = 0;
- }
- pop(b);
- errexit();
- lena = (bignum_length(a)>>2)-2;
- push(a);
- if ((int32)bignum_digits(b)[lenb] < 0)
- { b = negateb(b);
- /*
- * I MUST NOT do lenb = (bignum_length(b)>>2)-2; here because the
- * negateb may have failed and therefore handed back junk rather than
- * a bignum. E.g. an interrupt provoked by the user or a store-jam might
- * lead to a failure report here.
- */
- sign ^= 1;
- }
- else b = copyb(b);
- pop2(a, olda); /* original a, with original sign */
- errexit();
- lenb = (bignum_length(b)>>2)-2;
- /*
- * Now that the numbers are unsigned it could be that I can drop
- * a leading zero that had previously been necessary. That could reveal
- * that I have a one-word divisor after all...
- */
- if (bignum_digits(b)[lenb] == 0) lenb--;
- if (lenb == 0)
- { a = quotbn(a, bignum_digits(b)[0]);
- errexit();
- if ((sign & 2) != 0) nwork = -nwork;
- if ((sign & 1) != 0) a = negate(a);
- errexit();
- r = (Lisp_Object)(nwork & fix_mask);
- if (r == 0 || (int32)r == fix_mask) mv_2 = fixnum_of_int(nwork);
- else
- { push(a);
- if (signed_overflow(nwork))
- { if ((int32)nwork < 0)
- b = make_two_word_bignum(-1, clear_top_bit(nwork));
- else b = make_two_word_bignum(0, nwork);
- }
- else b = make_one_word_bignum(nwork);
- pop(a);
- errexit();
- mv_2 = b;
- }
- return a;
- }
- /*
- * Now the divisor is at least two digits long.
- */
- if (bignum_digits(a)[lena] == 0) lena--;
- /*
- * I will take a special case when the lengths of the two numbers
- * match since in that case the quotient will be only one digit long.
- * Also after I have filtered the lena<=lenb cases out, and have treated
- * one-word cases of b specially I will know that the numerator a has
- * at least two digits.
- */
- if (lena < lenb)
- { mv_2 = olda;
- return fixnum_of_int(0);
- }
- else if (lenb == lena)
- { unsigned32 p0 = bignum_digits(a)[lena-1], p1 = bignum_digits(a)[lena],
- q0 = bignum_digits(b)[lenb-1], q1 = bignum_digits(b)[lena],
- r0, r1, q, w, carry1;
- /*
- * If the dividend is smaller than the divisor I can return zero right now.
- */
- if (p1 < q1 || (p1 == q1 && p0 < q0))
- { mv_2 = olda;
- return fixnum_of_int(0);
- }
- /*
- * I scale the divisor so that it will have its top bit set (top wrt the
- * 31 bit field that is in use, that is), and scale the top two digits
- * of the dividend to match. The resulting values can then be divided to get
- * a pretty good estimate for the true quotient. Note that the division on
- * the next line is UNSIGNED arithmetic.
- */
- scale = 0x80000000U / ((unsigned32)1 + q1);
- Dmultiply(carry, w, q0, scale, 0);
- Dmultiply(q, w, q1, scale, carry);
- q = w;
- Dmultiply(carry, w, p0, scale, 0);
- Dmultiply(carry, w, p1, scale, carry);
- if (carry == q) q = 0x7fffffff;
- else
- { Ddivide(q, w, carry, w, q);
- q = w;
- }
- /*
- * q is now my estimated quotient, based on a quick look at high digits.
- */
- Dmultiply(carry, w, q0, q, 0);
- r0 = w;
- Dmultiply(carry, w, q1, q, carry);
- r1 = w;
- r0 = r0 - p0;
- if ((int32)r0 < 0)
- { r0 = clear_top_bit(r0);
- r1 = r1 - p1 - 1;
- }
- else r1 = r1 - p1;
- /* /*
- * the next line is a cop-out for now - if my estimated quotient
- * was close enough to the true value than the residual I get here
- * ought to be fairly small - if it is not I have bungled. Over a year
- * or so of tests I have not seen the disaster message triggered, but the
- * code should stay in until I can write the paragraph of comment that
- * should go here explaing why all is well...
- */
- if (carry != 0 && (int32)r1 < 0 &&
- carry != 1 && r1 != ~0x7fffffff)
- {
- err_printf("\n+++!!!+++ carry needed (line %d of \"%s\")\n",
- __LINE__, __FILE__);
- my_exit(EXIT_FAILURE);
- }
- /*
- * That obtained the remainder from (p1,p0)/(q1,q0) - now adjust q until
- * the remainder has the sign I want it to have.
- */
- /* /* I do not look at carry here - I have a nasty suspicion that I should.. */
- while ((int32)r1 > 0 ||
- (r1 == 0 && r0 != 0))
- { q--;
- r0 = r0 - q0;
- if ((int32)r0 < 0)
- { r0 = clear_top_bit(r0);
- r1 = r1 - q1 - 1;
- }
- else r1 = r1 - q1;
- }
- /*
- * Now q is (p1,p0)/(q1,q0), which is an over-estimate for the true
- * quotient, but by at worst 1. Compute a proper full-length remainder
- * now, using the original unscaled input numbers.
- */
- carry = 0;
- carry1 = 0xffffffffU;
- for (i=0; i<=lena; i++)
- { Dmultiply(carry, w, bignum_digits(b)[i], q, carry);
- carry1 = bignum_digits(a)[i] + clear_top_bit(~w) + top_bit(carry1);
- bignum_digits(a)[i] = clear_top_bit(carry1);
- }
- if ((int32)carry1 >= 0) /* need to correct it! */
- { q--;
- for (i=0; i<=lena; i++)
- { carry1 = bignum_digits(a)[i] + bignum_digits(b)[i] +
- top_bit(carry1);
- bignum_digits(a)[i] = clear_top_bit(carry1);
- }
- }
- /*
- * Now the quotient is correct - but since I want to hand back a neat
- * remainder I have more to do - I must convert the value stored in a
- * into a legitimate result, and store it in mv_2.
- */
- while (lena > 0 && bignum_digits(a)[lena] == 0) lena--;
- if ((bignum_digits(a)[lena] & 0x40000000) != 0) lena++;
- if (lena == 0) /* Maybe fixnum remainder? */
- { int32 rr = bignum_digits(a)[0];
- if (rr != 0 && (sign & 2) != 0)
- { rr = -rr;
- if ((rr & fix_mask) == fix_mask)
- { mv_2 = fixnum_of_int(rr);
- goto return_q;
- }
- }
- else if ((rr & fix_mask) == 0)
- { mv_2 = fixnum_of_int(rr);
- goto return_q;
- }
- }
- if ((sign & 2) != 0)
- { carry = 0xffffffffU;
- for (i=0; i<lena; i++)
- { carry = clear_top_bit(~bignum_digits(a)[i]) + top_bit(carry);
- bignum_digits(a)[i] = clear_top_bit(carry);
- }
- carry = ~bignum_digits(a)[i] + top_bit(carry);
- bignum_digits(a)[i] = carry;
- /*
- * if my remainder is -2^(31*n-1) then I can need to shrink lena here, which
- * seems like a messy case to have to consider.
- */
- if (carry == -1 && (bignum_digits(a)[i-1] & 0x40000000) != 0)
- { bignum_digits(a)[lena] = 0;
- lena--;
- bignum_digits(a)[i-1] |= ~0x7fffffff;
- }
- }
- /*
- * Now tidy up the space that I am about to discard, so that it will not
- * upset the garbage collector.
- */
- lenx = (bignum_length(a)>>2)-1;
- lenx |= 1; /* round up to allow for allocation in doublewords */
- lenc = (lena+1) | 1;
- if (lenc != lenx) /* space to discard? */
- bignum_digits(a)[lenc] = make_bighdr(lenx-lenc);
- numhdr(a) = make_bighdr(lena+2);
- mv_2 = a;
- return_q:
- if (q != 0 && (sign & 1) != 0)
- { q = -(int32)q;
- if ((q & fix_mask) == fix_mask) return fixnum_of_int(q);
- if ((q & 0x40000000) == 0)
- return make_two_word_bignum(-1, clear_top_bit(q));
- }
- else if ((q & fix_mask) == 0) return fixnum_of_int(q);
- else if ((q & 0x40000000) != 0) return make_two_word_bignum(0, q);
- return make_one_word_bignum(q);
- }
- /*
- * Now both a and b have been made positive, and their original signs
- * have been recorded so that I can tidy up at the end. The 1 bit in sign
- * tells me what sign the quotient will have, the 2 bit gives the sign
- * for the remainder. Both a and b have been copied so it will be OK to
- * overwrite them in the process of doing the division. The length of
- * a is strictly greater than that of b, which itself is at least 2 digits,
- * so this is a real case of long division, although the quotient may still
- * end up small (but non-zero).
- */
- /*
- * I find a scale factor to multiply a and b by so that the leading digit of
- * the divisor is fairly large. Again note that the arithmetic is UNSIGNED.
- */
- scale = 0x80000000U / ((unsigned32)1 + (unsigned32)bignum_digits(b)[lenb]);
- carry = 0;
- for (i=0; i<=lenb; i++)
- Dmultiply(carry, bignum_digits(b)[i],
- bignum_digits(b)[i], scale, carry);
- v1 = bignum_digits(b)[lenb];
- carry = 0;
- for (i=0; i<=lena; i++)
- Dmultiply(carry, bignum_digits(a)[i],
- bignum_digits(a)[i], scale, carry);
- /*
- * The scaled value of the numerator a may involve an extra digit
- * beyond those that I have space for in the vector. Hold this digit
- * in the variable 'carry'. Note that scaling the divisor did NOT cause
- * it to change length.
- */
- while (lena >= lenb)
- { unsigned32 w, h0, h1, v2, u2, carry1, carry2, qq, rr;
- if (carry == v1)
- { qq = 0x7fffffff;
- rr = carry + bignum_digits(a)[lena];
- }
- else
- { Ddivide(rr, w, carry, bignum_digits(a)[lena], v1);
- qq = w;
- }
- v2 = bignum_digits(b)[lenb-1];
- Dmultiply(h1, h0, v2, qq, 0);
- u2 = bignum_digits(a)[lena-1];
- if (h1 > rr ||
- (h1 == rr && h0 > u2))
- { qq--;
- h0 -= v2;
- if ((int32)h0 < 0)
- { h0 = clear_top_bit(h0);
- h1--;
- }
- rr += v1;
- if (h1 > rr ||
- (h1 == rr && h0 > u2)) qq--;
- }
- /*
- * Now subtract off q*b from a, up as far as lena.
- */
- carry1 = 0;
- carry2 = 0xffffffffU;
- for (i=0; i<=lenb; i++)
- { Dmultiply(carry1, w, bignum_digits(b)[i], qq, carry1);
- carry2 = bignum_digits(a)[lena-lenb+i] +
- clear_top_bit(~w) + top_bit(carry2);
- bignum_digits(a)[lena-lenb+i] = clear_top_bit(carry2);
- }
- carry2 = carry + clear_top_bit(~carry1) + top_bit(carry2);
- if ((int32)carry2 >= 0) /* overshot.. */
- { qq--;
- carry2 = 0;
- for (i=0; i<=lenb; i++)
- { carry2 = bignum_digits(a)[lena-lenb+i] +
- bignum_digits(b)[i] + top_bit(carry2);
- bignum_digits(a)[lena-lenb+i] = clear_top_bit(carry2);
- }
- }
- carry = bignum_digits(a)[lena];
- bignum_digits(a)[lena] = qq; /* good place to store the quotient */
- lena--;
- }
- lena = (bignum_length(a)>>2)-2;
- /*
- * Now the quotient is in the top digits of a, and the remainder is left
- * (scaled up) in the low bit (plus carry)
- */
- bignum_digits(b)[lenb] = carry/scale;
- carry = carry%scale;
- for (i=lenb-1; i>=0; i--)
- Ddivide(carry, bignum_digits(b)[i],
- carry, bignum_digits(a)[i], scale);
- if (carry != 0)
- { err_printf("+++!!!+++ Incorrect unscaling line %d file \"%s\" (%ld)\n",
- __LINE__, __FILE__, (long)carry);
- my_exit(EXIT_FAILURE);
- }
- for (i=0; i<=lena-lenb; i++)
- bignum_digits(a)[i] = bignum_digits(a)[i+lenb];
- for (;i < lena; i++) bignum_digits(a)[i] = 0;
- lena = lena - lenb;
- while (bignum_digits(a)[lena] == 0) lena--;
- /*
- * the above loop terminates properly since the quotient will be non-zero,
- * which in turn is because I was dividing p by q where p had strictly
- * more digits than q, and was hence strictly greater than q.
- */
- if ((bignum_digits(a)[lena] & 0x40000000) != 0) lena++;
- if ((sign & 1) != 0) /* quotient will be negative */
- { carry = 0xffffffffU;
- for (i=0; i<lena; i++)
- { carry = clear_top_bit(~bignum_digits(a)[i]) + top_bit(carry);
- bignum_digits(a)[i] = clear_top_bit(carry);
- }
- carry = ~bignum_digits(a)[i] + top_bit(carry);
- if (carry == -1 && (bignum_digits(a)[i-1] & 0x40000000) != 0)
- { bignum_digits(a)[lena] = 0;
- lena--;
- bignum_digits(a)[i-1] |= ~0x7fffffff;
- }
- else bignum_digits(a)[i] = carry;
- }
- lenx = (bignum_length(a)>>2)-1;
- lenx |= 1; /* round up to allow for allocation in doublewords */
- lenc = (lena+1) | 1;
- if (lenc != lenx) /* space to discard? */
- bignum_digits(a)[lenc] = make_bighdr(lenx-lenc);
- numhdr(a) = make_bighdr(lena+2);
- if (lena == 0)
- { int32 d = bignum_digits(a)[0];
- int32 x = d & fix_mask;
- if (x == 0 || x == fix_mask) a = fixnum_of_int(d);
- }
- /*
- * Here I need to tidy up the discarded space that could otherwise
- * kill the poor old garbage collector..
- */
- while (lenb >= 0 && bignum_digits(b)[lenb] == 0) lenb--;
- if (lenb < 0) mv_2 = fixnum_of_int(0); /* hah - zero remainder! */
- else
- {
- if ((bignum_digits(b)[lenb] & 0x40000000) != 0) lenb++;
- if ((sign & 2) != 0) /* need to negate the remainder */
- { carry = 0xffffffffU;
- for (i=0; i<lenb; i++)
- { carry = clear_top_bit(~bignum_digits(b)[i]) + top_bit(carry);
- bignum_digits(b)[i] = clear_top_bit(carry);
- }
- carry = ~bignum_digits(b)[i] + top_bit(carry);
- if (carry == -1 && (bignum_digits(b)[i-1] & 0x40000000) != 0)
- { bignum_digits(b)[lenb] = 0;
- lenb--;
- bignum_digits(b)[i-1] |= ~0x7fffffff;
- }
- else bignum_digits(b)[i] = carry;
- }
- lenx = (bignum_length(b)>>2)-1;
- lenx |= 1; /* round up to allow for allocation in doublewords */
- lenc = (lenb+1) | 1;
- if (lenc != lenx) /* space to discard? */
- bignum_digits(b)[lenc] = make_bighdr(lenx-lenc);
- numhdr(b) = make_bighdr(lenb+2);
- mv_2 = b;
- if (lenb == 0)
- { int32 r1, rr;
- rr = bignum_digits(b)[0];
- r1 = rr & fix_mask;
- if (r1 == 0 || r1 == fix_mask) mv_2 = fixnum_of_int(rr);
- }
- }
- return a;
- }
- #ifdef __powerc
- /* If you have trouble compiling this just comment it out, please */
- #pragma options(global_optimizer)
- #endif
- #ifdef COMMON
- #define quotbr(a, b) quotir(a, b)
- #define quotbc(a, b) quotic(a, b)
- #endif
- #define quotbf(a, b) quotsf(a, b)
- #ifdef COMMON
- static Lisp_Object quotri(Lisp_Object a, Lisp_Object b)
- {
- Lisp_Object w, nil;
- if (b == fixnum_of_int(1)) return a;
- else if (b == fixnum_of_int(0))
- return aerror2("bad arg for quotient", a, b);
- push3(a, b, C_nil);
- #define g stack[0]
- #define b stack[-1]
- #define a stack[-2]
- g = gcd(b, numerator(a));
- nil = C_nil;
- if (exception_pending()) goto fail;
- w = minusp(b);
- nil = C_nil;
- if (exception_pending()) goto fail;
- if (w)
- { g = negate(g); /* ensure denominator is +ve */
- nil = C_nil;
- if (exception_pending()) goto fail;
- }
- b = quot2(b, g);
- nil = C_nil;
- if (exception_pending()) goto fail;
- g = quot2(numerator(a), g);
- nil = C_nil;
- if (exception_pending()) goto fail;
- a = times2(b, denominator(a));
- nil = C_nil;
- if (exception_pending()) goto fail;
- w = make_ratio(g, a);
- popv(3);
- return w;
- fail:
- popv(3);
- return nil;
- #undef a
- #undef b
- #undef g
- }
- #define quotrs(a, b) quotsb(a, b)
- #define quotrb(a, b) quotib(a, b)
- static Lisp_Object quotrr(Lisp_Object a, Lisp_Object b)
- {
- Lisp_Object w, nil;
- push5(numerator(a), denominator(a),
- denominator(b), numerator(b), /* NB switched order */
- C_nil);
- #define g stack[0]
- #define db stack[-1]
- #define nb stack[-2]
- #define da stack[-3]
- #define na stack[-4]
- g = gcd(na, db);
- nil = C_nil;
- if (exception_pending()) goto fail;
- w = minusp(db);
- nil = C_nil;
- if (exception_pending()) goto fail;
- if (w)
- { g = negate(g);
- nil = C_nil;
- if (exception_pending()) goto fail;
- }
- na = quot2(na, g);
- nil = C_nil;
- if (exception_pending()) goto fail;
- db = quot2(db, g);
- nil = C_nil;
- if (exception_pending()) goto fail;
- g = gcd(nb, da);
- if (exception_pending()) goto fail;
- nil = C_nil;
- nb = quot2(nb, g);
- if (exception_pending()) goto fail;
- da = quot2(da, g);
- nil = C_nil;
- if (exception_pending()) goto fail;
- na = times2(na, nb);
- nil = C_nil;
- if (exception_pending()) goto fail;
- da = times2(da, db);
- nil = C_nil;
- if (exception_pending()) goto fail;
- w = make_ratio(na, da);
- popv(5);
- return w;
- fail:
- popv(5);
- return nil;
- #undef g
- #undef db
- #undef nb
- #undef da
- #undef na
- }
- #define quotrc(a, b) quotic(a, b)
- #define quotrf(a, b) quotsf(a, b)
- static Lisp_Object quotci(Lisp_Object a, Lisp_Object b)
- {
- Lisp_Object r = real_part(a), i = imag_part(a), nil;
- if (b == fixnum_of_int(0)) return aerror2("bad arg for quotient", a, b);
- push2(b, r);
- i = quot2(i, b);
- pop2(r, b);
- errexit();
- push(i);
- r = quot2(r, b);
- pop(i);
- errexit();
- return make_complex(r, i);
- }
- #define quotcs(a, b) quotci(a, b)
- #define quotcb(a, b) quotci(a, b)
- #define quotcr(a, b) quotci(a, b)
- #define quotcc(a, b) quotic(a, b)
- #define quotcf(a, b) quotci(a, b)
- #endif
- static Lisp_Object quotfi(Lisp_Object a, Lisp_Object b)
- {
- double d;
- if (b == fixnum_of_int(0)) return aerror2("bad arg for quotient", a, b);
- d = float_of_number(a) / (double)int_of_fixnum(b);
- return make_boxfloat(d, type_of_header(flthdr(a)));
- }
- static Lisp_Object quotfs(Lisp_Object a, Lisp_Object b)
- {
- double d = float_of_number(b);
- if (d == 0.0) return aerror2("bad arg for quotient", a, b);
- d = float_of_number(a) / d;
- return make_boxfloat(d, type_of_header(flthdr(a)));
- }
- #define quotfb(a, b) quotfs(a, b)
- #ifdef COMMON
- #define quotfr(a, b) quotfs(a, b)
- #define quotfc(a, b) quotic(a, b)
- #endif
- static Lisp_Object quotff(Lisp_Object a, Lisp_Object b)
- {
- #ifdef COMMON
- int32 ha = type_of_header(flthdr(a)), hb = type_of_header(flthdr(b));
- #endif
- double d = float_of_number(b);
- if (d == 0.0) return aerror2("bad arg for quotient", a, b);
- d = float_of_number(a) / d;
- #ifdef COMMON
- if (ha == TYPE_LONG_FLOAT || hb == TYPE_LONG_FLOAT)
- ha = TYPE_LONG_FLOAT;
- else if (ha == TYPE_DOUBLE_FLOAT || hb == TYPE_DOUBLE_FLOAT)
- ha = TYPE_DOUBLE_FLOAT;
- else ha = TYPE_SINGLE_FLOAT;
- return make_boxfloat(d, ha);
- #else
- return make_boxfloat(d, TYPE_DOUBLE_FLOAT);
- #endif
- }
- Lisp_Object quot2(Lisp_Object a, Lisp_Object b)
- {
- switch ((int)a & TAG_BITS)
- {
- case TAG_FIXNUM:
- switch ((int)b & TAG_BITS)
- {
- case TAG_FIXNUM:
- /*
- * This is where fixnum / fixnum arithmetic happens - the case I most want to
- * make efficient.
- */
- if (b == fixnum_of_int(0))
- return aerror2("bad arg for quotient", a, b);
- else
- { int32 r, aa, bb;
- aa = int_of_fixnum(a);
- bb = int_of_fixnum(b);
- /* calculate remainder and force its sign to be correct */
- r = aa % bb;
- if (aa < 0)
- { if (r > 0) r -= bb;
- }
- else if (r < 0) r += bb;
- /* then the division can be an exact one, as here */
- r = (aa - r)/bb;
- /*
- * the only case where dividing one fixnum by another can lead to
- * a bignum result is (-0x08000000/(-1)) which JUST overflows.
- */
- if (r != 0x08000000) return fixnum_of_int(r);
- else return make_one_word_bignum(r);
- }
- #ifdef COMMON
- case TAG_SFLOAT:
- return quotis(a, b);
- #endif
- case TAG_NUMBERS:
- { int32 hb = type_of_header(numhdr(b));
- switch (hb)
- {
- case TYPE_BIGNUM:
- return quotib(a, b);
- #ifdef COMMON
- case TYPE_RATNUM:
- return quotir(a, b);
- case TYPE_COMPLEX_NUM:
- return quotic(a, b);
- #endif
- default:
- return aerror1("bad arg for quotient", b);
- }
- }
- case TAG_BOXFLOAT:
- return quotif(a, b);
- default:
- return aerror1("bad arg for quotient", b);
- }
- #ifdef COMMON
- case TAG_SFLOAT:
- switch ((int)b & TAG_BITS)
- {
- case TAG_FIXNUM:
- return quotsi(a, b);
- case TAG_SFLOAT:
- { Float_union aa, bb;
- aa.i = a - TAG_SFLOAT;
- bb.i = b - TAG_SFLOAT;
- aa.f = (float) (aa.f / bb.f);
- return (aa.i & ~(int32)0xf) + TAG_SFLOAT;
- }
- case TAG_NUMBERS:
- { int32 hb = type_of_header(numhdr(b));
- switch (hb)
- {
- case TYPE_BIGNUM:
- return quotsb(a, b);
- case TYPE_RATNUM:
- return quotsr(a, b);
- case TYPE_COMPLEX_NUM:
- return quotsc(a, b);
- default:
- return aerror1("bad arg for quotient", b);
- }
- }
- case TAG_BOXFLOAT:
- return quotsf(a, b);
- default:
- return aerror1("bad arg for quotient", b);
- }
- #endif
- case TAG_NUMBERS:
- { int32 ha = type_of_header(numhdr(a));
- switch (ha)
- {
- case TYPE_BIGNUM:
- switch ((int)b & TAG_BITS)
- {
- case TAG_FIXNUM:
- return quotbi(a, b);
- #ifdef COMMON
- case TAG_SFLOAT:
- return quotbs(a, b);
- #endif
- case TAG_NUMBERS:
- { int32 hb = type_of_header(numhdr(b));
- switch (hb)
- {
- case TYPE_BIGNUM:
- return quotbb(a, b);
- #ifdef COMMON
- case TYPE_RATNUM:
- return quotbr(a, b);
- case TYPE_COMPLEX_NUM:
- return quotbc(a, b);
- #endif
- default:
- return aerror1("bad arg for quotient", b);
- }
- }
- case TAG_BOXFLOAT:
- return quotbf(a, b);
- default:
- return aerror1("bad arg for quotient", b);
- }
- #ifdef COMMON
- case TYPE_RATNUM:
- switch ((int)b & TAG_BITS)
- {
- case TAG_FIXNUM:
- return quotri(a, b);
- case TAG_SFLOAT:
- return quotrs(a, b);
- case TAG_NUMBERS:
- { int32 hb = type_of_header(numhdr(b));
- switch (hb)
- {
- case TYPE_BIGNUM:
- return quotrb(a, b);
- case TYPE_RATNUM:
- return quotrr(a, b);
- case TYPE_COMPLEX_NUM:
- return quotrc(a, b);
- default:
- return aerror1("bad arg for quotient", b);
- }
- }
- case TAG_BOXFLOAT:
- return quotrf(a, b);
- default:
- return aerror1("bad arg for quotient", b);
- }
- case TYPE_COMPLEX_NUM:
- switch ((int)b & TAG_BITS)
- {
- case TAG_FIXNUM:
- return quotci(a, b);
- case TAG_SFLOAT:
- return quotcs(a, b);
- case TAG_NUMBERS:
- { int32 hb = type_of_header(numhdr(b));
- switch (hb)
- {
- case TYPE_BIGNUM:
- return quotcb(a, b);
- case TYPE_RATNUM:
- return quotcr(a, b);
- case TYPE_COMPLEX_NUM:
- return quotcc(a, b);
- default:
- return aerror1("bad arg for quotient", b);
- }
- }
- case TAG_BOXFLOAT:
- return quotcf(a, b);
- default:
- return aerror1("bad arg for quotient", b);
- }
- #endif
- default: return aerror1("bad arg for quotient", a);
- }
- }
- case TAG_BOXFLOAT:
- switch ((int)b & TAG_BITS)
- {
- case TAG_FIXNUM:
- return quotfi(a, b);
- #ifdef COMMON
- case TAG_SFLOAT:
- return quotfs(a, b);
- #endif
- case TAG_NUMBERS:
- { int32 hb = type_of_header(numhdr(b));
- switch (hb)
- {
- case TYPE_BIGNUM:
- return quotfb(a, b);
- #ifdef COMMON
- case TYPE_RATNUM:
- return quotfr(a, b);
- case TYPE_COMPLEX_NUM:
- return quotfc(a, b);
- #endif
- default:
- return aerror1("bad arg for quotient", b);
- }
- }
- case TAG_BOXFLOAT:
- return quotff(a, b);
- default:
- return aerror1("bad arg for quotient", b);
- }
- default:
- return aerror1("bad arg for quotient", a);
- }
- }
- #ifdef COMMON
- /*
- * The above version of quot2 is as required for the Standard Lisp QUOTIENT
- * function, and it is called from various internal places in CSL/CCL, eg
- * from within the code for TRUNCATE. Next I have something that will be very
- * similar in behaviour, but which turns quotients of integers into
- * rational numbers when that is necessary.
- */
- /***** Not reconstructed yet!! */
- Lisp_Object CLquot2(Lisp_Object a, Lisp_Object b)
- {
- switch ((int)a & TAG_BITS)
- {
- case TAG_FIXNUM:
- switch ((int)b & TAG_BITS)
- {
- case TAG_FIXNUM:
- /*
- * This is where fixnum / fixnum arithmetic happens - the case I most want to
- * make efficient.
- */
- if (b == fixnum_of_int(0))
- return aerror2("bad arg for /", a, b);
- else
- { int32 r, aa, bb, w;
- aa = int_of_fixnum(a);
- bb = int_of_fixnum(b);
- if (bb < 0) aa = -aa, bb = -bb;
- r = aa % bb;
- if (r == 0) /* Exact division case */
- { r = aa/bb;
- /*
- * the only case where dividing one fixnum by another can lead to
- * a bignum result is (-0x08000000/(-1)) which JUST overflows.
- */
- if (r != 0x08000000) return fixnum_of_int(r);
- else return make_one_word_bignum(r);
- }
- w = bb;
- if (r < 0) r = -r;
- while (r != 0)
- { int32 w1 = w % r;
- w = r;
- r = w1;
- }
- aa = aa / w;
- bb = bb / w;
- return make_ratio(fixnum_of_int(aa), fixnum_of_int(bb));
- }
- #ifdef COMMON
- case TAG_SFLOAT:
- return quotis(a, b);
- #endif
- case TAG_NUMBERS:
- { int32 hb = type_of_header(numhdr(b));
- switch (hb)
- {
- case TYPE_BIGNUM:
- return CLquotib(a, b);
- #ifdef COMMON
- case TYPE_RATNUM:
- return quotir(a, b);
- case TYPE_COMPLEX_NUM:
- return quotic(a, b);
- #endif
- default:
- return aerror1("bad arg for /", b);
- }
- }
- case TAG_BOXFLOAT:
- return quotif(a, b);
- default:
- return aerror1("bad arg for /", b);
- }
- #ifdef COMMON
- case TAG_SFLOAT:
- switch ((int)b & TAG_BITS)
- {
- case TAG_FIXNUM:
- return quotsi(a, b);
- case TAG_SFLOAT:
- { Float_union aa, bb;
- aa.i = a - TAG_SFLOAT;
- bb.i = b - TAG_SFLOAT;
- aa.f = (float) (aa.f / bb.f);
- return (aa.i & ~(int32)0xf) + TAG_SFLOAT;
- }
- case TAG_NUMBERS:
- { int32 hb = type_of_header(numhdr(b));
- switch (hb)
- {
- case TYPE_BIGNUM:
- return quotsb(a, b);
- case TYPE_RATNUM:
- return quotsr(a, b);
- case TYPE_COMPLEX_NUM:
- return quotsc(a, b);
- default:
- return aerror1("bad arg for /", b);
- }
- }
- case TAG_BOXFLOAT:
- return quotsf(a, b);
- default:
- return aerror1("bad arg for /", b);
- }
- #endif
- case TAG_NUMBERS:
- { int32 ha = type_of_header(numhdr(a));
- switch (ha)
- {
- case TYPE_BIGNUM:
- switch ((int)b & TAG_BITS)
- {
- case TAG_FIXNUM:
- return CLquotbi(a, b);
- #ifdef COMMON
- case TAG_SFLOAT:
- return quotbs(a, b);
- #endif
- case TAG_NUMBERS:
- { int32 hb = type_of_header(numhdr(b));
- switch (hb)
- {
- case TYPE_BIGNUM:
- return CLquotbb(a, b);
- #ifdef COMMON
- case TYPE_RATNUM:
- return quotbr(a, b);
- case TYPE_COMPLEX_NUM:
- return quotbc(a, b);
- #endif
- default:
- return aerror1("bad arg for /", b);
- }
- }
- case TAG_BOXFLOAT:
- return quotbf(a, b);
- default:
- return aerror1("bad arg for /", b);
- }
- #ifdef COMMON
- case TYPE_RATNUM:
- switch ((int)b & TAG_BITS)
- {
- case TAG_FIXNUM:
- return quotri(a, b);
- case TAG_SFLOAT:
- return quotrs(a, b);
- case TAG_NUMBERS:
- { int32 hb = type_of_header(numhdr(b));
- switch (hb)
- {
- case TYPE_BIGNUM:
- return quotrb(a, b);
- case TYPE_RATNUM:
- return quotrr(a, b);
- case TYPE_COMPLEX_NUM:
- return quotrc(a, b);
- default:
- return aerror1("bad arg for /", b);
- }
- }
- case TAG_BOXFLOAT:
- return quotrf(a, b);
- default:
- return aerror1("bad arg for /", b);
- }
- case TYPE_COMPLEX_NUM:
- switch ((int)b & TAG_BITS)
- {
- case TAG_FIXNUM:
- return quotci(a, b);
- case TAG_SFLOAT:
- return quotcs(a, b);
- case TAG_NUMBERS:
- { int32 hb = type_of_header(numhdr(b));
- switch (hb)
- {
- case TYPE_BIGNUM:
- return quotcb(a, b);
- case TYPE_RATNUM:
- return quotcr(a, b);
- case TYPE_COMPLEX_NUM:
- return quotcc(a, b);
- default:
- return aerror1("bad arg for /", b);
- }
- }
- case TAG_BOXFLOAT:
- return quotcf(a, b);
- default:
- return aerror1("bad arg for /", b);
- }
- #endif
- default: return aerror1("bad arg for /", a);
- }
- }
- case TAG_BOXFLOAT:
- switch ((int)b & TAG_BITS)
- {
- case TAG_FIXNUM:
- return quotfi(a, b);
- #ifdef COMMON
- case TAG_SFLOAT:
- return quotfs(a, b);
- #endif
- case TAG_NUMBERS:
- { int32 hb = type_of_header(numhdr(b));
- switch (hb)
- {
- case TYPE_BIGNUM:
- return quotfb(a, b);
- #ifdef COMMON
- case TYPE_RATNUM:
- return quotfr(a, b);
- case TYPE_COMPLEX_NUM:
- return quotfc(a, b);
- #endif
- default:
- return aerror1("bad arg for /", b);
- }
- }
- case TAG_BOXFLOAT:
- return quotff(a, b);
- default:
- return aerror1("bad arg for /", b);
- }
- default:
- return aerror1("bad arg for /", a);
- }
- }
- #endif
- /* end of arith03.c */
|