/* arith02.c Copyright (C) 1990-2002 Codemist Ltd */
/*
* Arithmetic functions.
* Multiplication of generic numbers
* and in particular a lot of bignum support.
*
*/
/*
* This code may be used and modified, and redistributed in binary
* or source form, subject to the "CCL Public License", which should
* accompany it. This license is a variant on the BSD license, and thus
* permits use of code derived from this in either open and commercial
* projects: but it does require that updates to this code be made
* available back to the originators of the package.
* Before merging other code in with this or linking this code
* with other packages or libraries please check that the license terms
* of the other material are compatible with those of this.
*/
/* Signature: 1884c4f0 08-Apr-2002 */
#include <stdarg.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#ifdef __WATCOMC__
#include <float.h>
#endif
#include "machine.h"
#include "tags.h"
#include "cslerror.h"
#include "externs.h"
#include "arith.h"
#ifdef TIMEOUT
#include "timeout.h"
#endif
/*
* Now for multiplication
*/
/*
* I provide symbols IMULTIPLY and IDIVIDE which can be asserted if the
* corresponding routines have been provided elsewhere (e.g. in machine
* code for extra speed)
*/
#ifndef IMULTIPLY
#ifdef MULDIV64
/*
* If MULDIV64 is asserted this function in not in fact needed
* since a macro in arith.h arranges that the multiplication is done
* in-line. However this version is left here in the source code as
* convenient documentation of what the function needs to do.
*/
unsigned32 Imultiply(unsigned32 *rlow, unsigned32 a, unsigned32 b, unsigned32 c)
/*
* (result, *rlow) = a*b + c as 62-bit product
*
* rlow may well be the same as one of a, b or c so the
* assignment to *rlow must not be done until the calculation is
* complete. Inputs a and b are 31-bit values (i.e. they have
* their top bit zero and are to be treated as positive value.
* c may use all 32 bits, but again is treated as unsigned.
* The result is computed and the low 31 bits placed in rlow,
* (the top bit of rlow will end up zero) and the top part
* of the result is returned.
*/
{
/* NB the value r cound be int64 or unsigned64 - it does not matter */
unsigned64 r = (unsigned64)a*(unsigned64)b + (unsigned64)c;
*rlow = (unsigned32)(r & 0x7fffffff);
return (unsigned32)(r >> 31);
}
#else
#ifdef OLDER_VERSION
unsigned32 Imultiply(unsigned32 *rlow, unsigned32 a, unsigned32 b, unsigned32 c)
/*
* (result, *rlow) = a*b + c as 62-bit product
*
* The code given here forms the produce by doing four single-precision
* (16*16->32) multiplies and using shifts etc to glue the partial results
* together.
*/
{
unsigned32 ah = a >> 16, bh = b >> 16, ch = c >> 16;
unsigned32 w0, w1, w2, w3, w4;
a -= ah << 16;
b -= bh << 16;
c -= ch << 16; /* a, b and c now split into high/low parts */
w0 = a*b; /* One... */
w1 = w0 >> 16;
w0 = w0 - (w1 << 16) + c;
w1 += ch;
w2 = a*bh; /* Two... */
w3 = w2 >> 16;
w1 += w2 - (w3 << 16);
w2 = ah*b; /* Three... */
w4 = w2 >> 16;
w1 += w2 - (w4 << 16);
w2 = w0 >> 16;
w1 += w2;
w0 -= w2 << 16;
w2 = ah*bh + w3 + w4; /* Four 16-bit multiplies done in all. */
w2 += w1 >> 16;
/* Here I do a minor shift to split the result at bit 31 rather than 32 */
*rlow = w0 + ((w1 & 0x7fff) << 16);
return (w2<<1) + ((w1>>15) & 1);
}
#else
unsigned32 Imultiply(unsigned32 *rlow, unsigned32 a, unsigned32 b, unsigned32 c)
/*
* (result, *rlow) = a*b + c as 62-bit product
*
* The code given here forms the produce by doing three single-precision
* (16*16->32) multiplies and using shifts etc to glue the partial results
* together. This is slightly faster than the above (maybe simpler?) code
* on at least some machines.
*/
{
unsigned32 ah, bh;
unsigned32 w0, w1, w2;
ah = a >> 16;
/*
* On some machines I know that multi-bit shifts are especially painful,
* while on others it is nasty to access the literal value needed as a
* mask here. Hence I make some show of providing an alternative.
*/
#ifdef FAST_SHIFTS
a -= ah << 16;
#else
a &= 0xffff;
#endif
bh = b >> 16;
#ifdef FAST_SHIFTS
b -= bh << 16;
#else
b &= 0xffff;
#endif
/*
* At present I can not see any way of issuing any of these multiplies
* any earlier, or of doing useful work between the times that I launch
* them, or even delaying before I use their results. This will be rather
* sad on machines where multiplication could overlap with other operations.
*/
w2 = ah*bh;
w1 = (a + ah)*(b + bh);
w0 = a*b;
/*
* The largest exact result that can be computed by the next line (given
* that a and b start off just 31 bit) is 0xfffd0002
*/
w1 = (w1 - w2) - w0; /* == (a*bh + b*ah) */
/*
* I split into 30 bits in the lower word and 32 in the upper, so I have
* 2 bits available for temporary carry effects in the lower part.
*/
w2 = (w2 << 2) + (w1 >> 14) + (w0 >> 30) + (c >> 30);
w1 &= 0x3fff;
w0 = (c & 0x3fffffff) + (w0 & 0x3fffffffU) + (w1 << 16);
w0 += ((w2 & 1) << 30);
w2 = (w2 >> 1) + (w0 >> 31);
*rlow = w0 & 0x7fffffffU;
return w2;
}
#endif
#endif
#endif /* IMULTIPLY */
static Lisp_Object timesii(Lisp_Object a, Lisp_Object b)
/*
* multiplying two fixnums together is much messier than adding them,
* mainly because the result can easily be a two-word bignum
*/
{
unsigned32 aa = (unsigned32)int_of_fixnum(a),
bb = (unsigned32)int_of_fixnum(b);
unsigned32 temp, low, high;
/*
* Multiplication by 0 or by -1 is just possibly common enough to be worth
* filtering out the following special cases. Avoidance of the tedious
* checks for overflow may make this useful even if Imultiply is very fast.
*/
if (aa <= 1)
{ if (aa == 0) return fixnum_of_int(0);
else return b;
}
else if (bb <= 1)
{ if (bb == 0) return fixnum_of_int(0);
else return a;
}
/*
* I dump the low part of the product in temp then copy to a variable low
* because temp has to have its address taken and so is not a candidate for
* living in a register.
*/
Dmultiply(high, temp, clear_top_bit(aa), clear_top_bit(bb), 0);
/*
* The result handed back here has only 31 bits active in its low part.
*/
low = temp;
/*
* The next two lines convert the unsigned product produced by Imultiply
* into a signed product.
*/
if ((int32)aa < 0) high -= bb;
if ((int32)bb < 0) high -= aa;
if ((high & 0x40000000) == 0) high = clear_top_bit(high);
else high = set_top_bit(high);
if (high == 0 && (low & 0x40000000) == 0)
{ /* one word positive result */
if (low <= 0x07ffffff) return fixnum_of_int(low);
else return make_one_word_bignum(low);
}
else if (high == -1 && (low & 0x40000000) != 0)
{ /* one word negative result */
low = set_top_bit(low);
if ((low & fix_mask) == fix_mask) return fixnum_of_int(low);
else return make_one_word_bignum(low);
}
else return make_two_word_bignum(high, low);
}
#ifdef COMMON
static Lisp_Object timesis(Lisp_Object a, Lisp_Object b)
{
Float_union bb;
bb.i = b - TAG_SFLOAT;
bb.f = (float) ((float)int_of_fixnum(a) * bb.f);
return (bb.i & ~(int32)0xf) + TAG_SFLOAT;
}
#endif
static Lisp_Object timesib(Lisp_Object a, Lisp_Object b)
{
int32 aa = int_of_fixnum(a), lenb, i;
unsigned32 carry, ms_dig, w;
Lisp_Object c, nil;
/*
* I will split off the (easy) cases of the fixnum being -1, 0 or 1.
*/
if (aa == 0) return fixnum_of_int(0);
else if (aa == 1) return b;
else if (aa == -1) return negateb(b);
lenb = bignum_length(b);
push(b);
c = getvector(TAG_NUMBERS, TYPE_BIGNUM, lenb);
pop(b);
errexit();
lenb = (lenb-CELL)/4;
if (aa < 0)
{ aa = -aa;
carry = 0xffffffffU;
for (i=0; i<lenb-1; i++)
{ carry = clear_top_bit(~bignum_digits(b)[i]) + top_bit(carry);
bignum_digits(c)[i] = clear_top_bit(carry);
}
/*
* I do the most significant digit separately.
*/
carry = clear_top_bit(~bignum_digits(b)[i]) + top_bit(carry);
/*
* there is a special case needed here - if b started off as a number
* like 0xc0000000,0,0,0 then negating it would call for extension of
* the bignum c (this is the usual assymetry in the range of twos-
* complement numbers). But if I detect that case specially I can
* observe that the ONLY case where negation overflows is where the
* negated value is exactly a power of 2 such .. an easy thing to
* multiply a by! Furthermore the power of two involved is known to
* by 30 mod 31.
*/
if (carry == 0x40000000)
{ bignum_digits(c)[i] = (aa & 1) << 30;
aa = aa >> 1;
goto extend_by_one_word;
}
if ((carry & 0x40000000) != 0) carry = set_top_bit(carry);
bignum_digits(c)[i] = carry;
}
else
{ for (i=0; i<lenb; i++) bignum_digits(c)[i] = bignum_digits(b)[i];
}
/*
* Now c is a copy of b (negated if necessary) and I just want to
* multiply it by the positive value a. This is the heart of the
* procedure. Re-write Imultiply in assembly code if you want it
* to go faster. See the top of this file for a portable sample
* implementation of Imultiply, which gives back a result as a pair
* of 31-bit values. But note that if the computer has been marked as
* supporting a 64-bit integer data-type then Imultiply macro-expands
* into code that uses 64-bit intermediate values.
*/
carry = 0;
/*
* here aa is > 0, and a fortiori the 0x80000000 bit of aa is clear,
* so I do not have to worry about the difference between 31 and 32
* bit values for aa.
*/
for (i=0; i<lenb-1; i++)
Dmultiply(carry, bignum_digits(c)[i], bignum_digits(c)[i],
(unsigned32)aa, carry);
ms_dig = bignum_digits(c)[i];
Dmultiply(carry, w, clear_top_bit(ms_dig), (unsigned32)aa, carry);
/*
* After forming the product to (lenb) digits I need to see if there
* is any overflow. Calculate what the next digit would be, sign
* extending into the 0x80000000 bit as necessary.
*/
if ((carry & 0x40000000) != 0) carry = set_top_bit(carry);
if (((int32)ms_dig) < 0) carry -= aa;
aa = (int32)carry;
if (aa == -1 && (w & 0x40000000) != 0)
{ bignum_digits(c)[i] = set_top_bit(w);
return c;
}
bignum_digits(c)[i] = w;
if (aa == 0 && (w & 0x40000000) == 0) return c;
/*
* drop through to extend the number by a word - note that because I
* am multiplying by a fixnum it is only possible to have to expand by
* just one word.
*/
extend_by_one_word:
#ifdef ADDRESS_64
if ((lenb & 1) != 0)
#else
if ((lenb & 1) == 0)
#endif
/*
* Here there was a padder word that I can expand into.
*/
{ bignum_digits(c)[lenb] = aa;
numhdr(c) += pack_hdrlength(1);
return c;
}
/*
* Need to allocate more space to grow into. I need to grow by just 4 bytes.
*/
push(c);
a = getvector(TAG_NUMBERS, TYPE_BIGNUM, CELL+4+4*lenb);
pop(c);
errexit();
for (i=0; i<lenb; i++)
bignum_digits(a)[i] = bignum_digits(c)[i];
bignum_digits(a)[i] = aa;
bignum_digits(a)[i+1] = 0; /* the padder word */
return a;
}
#ifdef COMMON
static Lisp_Object timesir(Lisp_Object a, Lisp_Object b)
/*
* multiply integer (fixnum or bignum) by ratio.
*/
{
Lisp_Object nil = C_nil;
Lisp_Object w = nil;
if (a == fixnum_of_int(0)) return a;
else if (a == fixnum_of_int(1)) return b;
push3(b, a, nil);
#define g stack[0]
#define a stack[-1]
#define b stack[-2]
g = gcd(a, denominator(b));
nil = C_nil;
if (exception_pending()) goto fail;
a = quot2(a, g);
nil = C_nil;
if (exception_pending()) goto fail;
g = quot2(denominator(b), g);
nil = C_nil;
if (exception_pending()) goto fail;
a = times2(a, numerator(b));
nil = C_nil;
if (exception_pending()) goto fail;
/*
* make_ratio tidies things up if the denominator was exactly 1
*/
w = make_ratio(a, g);
fail:
popv(3);
return w;
#undef a
#undef b
#undef g
}
static Lisp_Object timesic(Lisp_Object a, Lisp_Object b)
/*
* multiply an arbitrary non-complex number by a complex one
*/
{
Lisp_Object nil;
Lisp_Object r = real_part(b), i = imag_part(b);
push2(a, r);
i = times2(a, i);
pop2(r, a);
errexit();
push(i);
r = times2(a, r);
pop(i);
errexit();
return make_complex(r, i);
}
#endif
static Lisp_Object timesif(Lisp_Object a, Lisp_Object b)
{
double d = (double)int_of_fixnum(a) * float_of_number(b);
return make_boxfloat(d, type_of_header(flthdr(b)));
}
#ifdef COMMON
#define timessi(a, b) timesis(b, a)
static Lisp_Object timessb(Lisp_Object a, Lisp_Object b)
{
double d = float_of_number(a) * float_of_number(b);
return make_sfloat(d);
}
#define timessr(a, b) timessb(a, b)
#define timessc(a, b) timesic(a, b)
#endif
static Lisp_Object timessf(Lisp_Object a, Lisp_Object b)
{
double d = float_of_number(a) * float_of_number(b);
return make_boxfloat(d, type_of_header(flthdr(b)));
}
#define timesbi(a, b) timesib(b, a)
#ifdef COMMON
#define timesbs(a, b) timessb(b, a)
#endif
/*
* Now for bignum multiplication - made more than comfortably complicated
* by a desire to make it go fast for very big numbers.
*/
#ifndef KARATSUBA_CUTOFF
/*
* I have conducted some experiments on one machine to find out what the
* best cut-off value here is. The exact value chosen is not very
* critical, and the fancy techniques do not pay off until numbers get
* a lot bigger than this length (which is expressed in units of 31-bit
* words, i.e. about 10 decimals). Anyone who wants may recompile with
* alternative values to try to get the system fine-tuned for their
* own computer - but I do not expect it to be possible to achieve much
* by so doing.
*/
#define KARATSUBA_CUTOFF 12
#endif
static void long_times(unsigned32 *c, unsigned32 *a, unsigned32 *b,
unsigned32 *d, int32 lena, int32 lenb, int32 lenc);
static void long_times1(unsigned32 *c, unsigned32 *a, unsigned32 *b,
unsigned32 *d, int32 lena, int32 lenb, int32 lenc)
/*
* Here both a and b are big, with lena <= lenb. Split each into two chunks
* of size (lenc/4), say (a1,a2) and (b1,b2), and compute each of
* a1*b1
* a2*b2
* (a1+a2)*(b1+b2)
* where in the last case a1+a2 and b1+b2 can be computed keeping their
* carry bits as k1 and k2 (so that a1+a2 and b1+b2 are restricted to
* size lenb). The chunks can then be combined with a few additions and
* subtractions to form the product that is really wanted. If in fact a
* was shorter than lenc/4 (so that after a is split up the top half is
* all zero) I do things in a more straightforward way. I require that on
* entry to this code lenc<4 < lenb <= lenc/2.
*/
{
int32 h = lenc/4; /* lenc must have been made even enough... */
int32 lena1 = lena - h;
int32 lenb1 = lenb - h;
unsigned32 carrya, carryb;
int32 i;
/*
* if the top half of a would be all zero I go through a separate path,
* doing just two subsidiary multiplications.
*/
if (lena1 <= 0)
{ long_times(c+h, a, b+h, d, lena, lenb-h, 2*h);
for (i=0; i<h; i++) c[3*h+i] = 0;
long_times(d, a, b, c, lena, h, 2*h);
for (i=0; i<h; i++) c[i] = d[i];
carrya = 0;
for (;i<2*h; i++)
{ unsigned32 w = c[i] + d[i] + carrya;
c[i] = clear_top_bit(w);
carrya = w >> 31;
}
for (;carrya!=0;i++)
{ unsigned32 w = c[i] + 1;
c[i] = clear_top_bit(w);
carrya = w >> 31;
}
return;
}
/*
* form (a1+a2) and (b1+b2).
*/
carrya = 0;
for (i=0; i<h; i++)
{ unsigned32 w = a[i] + carrya;
if (i < lena1) w += a[h+i];
d[i] = clear_top_bit(w);
carrya = w >> 31;
}
carryb = 0;
for (i=0; i<h; i++)
{ unsigned32 w = b[i] + carryb;
if (i < lenb1) w += b[h+i];
d[h+i] = clear_top_bit(w);
carryb = w >> 31;
}
long_times(c+h, d, d+h, c, h, h, 2*h);
/*
* Adjust to allow for the cases of a1+a2 or b1+b2 overflowing
* by a single bit.
*/
c[3*h] = carrya & carryb;
if (carrya != 0)
{ carrya = 0;
for (i=0; i<h; i++)
{ unsigned32 w = c[2*h+i] + d[h+i] + carrya;
c[2*h+i] = clear_top_bit(w);
carrya = w >> 31;
}
}
if (carryb != 0)
{ carryb = 0;
for (i=0; i<h; i++)
{ unsigned32 w = c[2*h+i] + d[i] + carryb;
c[2*h+i] = clear_top_bit(w);
carryb = w >> 31;
}
}
c[3*h] += carrya + carryb;
for (i=1; i<h; i++) c[3*h+i] = 0;
/*
* Now (a1+a2)*(b1+b2) should have been computed totally properly
*/
for (i=0; i<h; i++) d[h+i] = 0;
/*
* multiply out a1*b1, where note that a1 and b1 may be less long
* than h, but not by much.
*/
long_times(d, a+h, b+h, c, lena-h, lenb-h, 2*h);
carrya = 0;
for (i=0; i<2*h; i++)
{ unsigned32 w = c[2*h+i] + d[i] + carrya;
c[2*h+i] = clear_top_bit(w);
carrya = w >> 31;
}
carrya = 0;
for (i=0; i<2*h; i++)
{ unsigned32 w = c[h+i] - d[i] - carrya;
c[h+i] = clear_top_bit(w);
carrya = w >> 31;
}
for (; carrya!=0 && i<3*h; i++)
{ unsigned32 w = c[h+i] - 1;
c[h+i] = clear_top_bit(w);
carrya = w >> 31;
}
/*
* multiply out a2*b2
*/
long_times(d, a, b, c, h, h, 2*h);
for (i=0; i<h; i++) c[i] = d[i];
carrya = 0;
for (; i<2*h; i++)
{ unsigned32 w = c[i] + d[i] + carrya;
c[i] = clear_top_bit(w);
carrya = w >> 31;
}
for (; carrya!=0 && i<4*h; i++)
{ unsigned32 w = c[i] + 1;
c[i] = clear_top_bit(w);
carrya = w >> 31;
}
carrya = 0;
for (i=0; i<2*h; i++)
{ unsigned32 w = c[h+i] - d[i] - carrya;
c[h+i] = clear_top_bit(w);
carrya = w >> 31;
}
for (; carrya!=0 && i<3*h; i++)
{ unsigned32 w = c[h+i] - 1;
c[h+i] = clear_top_bit(w);
carrya = w >> 31;
}
/*
* The product is now complete
*/
}
static void long_times2(unsigned32 *c, unsigned32 *a, unsigned32 *b,
int32 lena, int32 lenb, int32 lenc)
/*
* This case is standard old fashioned long multiplication. Dump the
* result into c.
*/
{
int32 i;
for (i=0; i<lenc; i++) c[i] = 0;
for (i=0; i<lena; i++)
{ unsigned32 carry = 0, da = a[i];
int32 j;
/*
* When I multiply by (for instance) a high power of 2 there will
* be plenty of zero digits in the number being worked with, and
* so the test da!=0 will save something useful.
*/
if (da != 0)
{ for (j=0; j<lenb; j++)
{ int32 k = i + j;
Dmultiply(carry, c[k], da, b[j],
/* NB the addition here is OK and fits into a 32-bit unsigned result */
carry + c[k]);
}
c[i+j] = carry;
}
}
}
static void long_times(unsigned32 *c, unsigned32 *a, unsigned32 *b,
unsigned32 *d, int32 lena, int32 lenb, int32 lenc)
/*
* This decides if a multiplication is big enough to benefit from
* decomposition a la Karatsuba.
* In recursive entries through here out of long_times1() the numbers a
* and b may have shrunk in ways that mean I need to reconsider the
* precision to which I am working. This must leave c filled out all
* the way to lenc, with padding 0s if necessary.
*/
{
if (lenb < lena)
{ unsigned32 *t1;
int32 t2;
t1 = a; a = b; b = t1;
t2 = lena; lena = lenb; lenb = t2;
}
if (4*lenb <= lenc) /* In this case I should shrink lenc a bit.. */
{ int32 newlenc = (lenb+1)/2;
int k = 0;
while (newlenc > KARATSUBA_CUTOFF)
{ newlenc = (newlenc + 1)/2;
k++;
}
while (k != 0)
{ newlenc = 2*newlenc;
k--;
}
newlenc = 4*newlenc;
while (lenc > newlenc) c[--lenc] = 0;
}
if (lena > KARATSUBA_CUTOFF) long_times1(c, a, b, d, lena, lenb, lenc);
else long_times2(c, a, b, lena, lenb, lenc);
}
static Lisp_Object timesbb(Lisp_Object a, Lisp_Object b)
/*
* a and b are both guaranteed to be bignums when I call this
* procedure.
*/
{
int sign = 1;
Lisp_Object c, d, nil;
int32 lena, lenb, lenc, i;
lena = (bignum_length(a) - CELL)/4;
lenb = (bignum_length(b) - CELL)/4;
if (lena == 1 && lenb == 1)
/*
* I am going to deem multiplication of two one-word bignums worthy of
* a special case, since it is probably fairly common and it will be cheap
* enough that avoiding overheads might matter. I still need to split
* off the signs, because Imultiply can only deal with 31-bit unsigned values.
* One-word bignums each have value at least 2^27 (or else they would have
* been represented as fixnums) so the result here will always be a two-word
* bignum.
*/
{ int32 va = (int32)bignum_digits(a)[0],
vb = (int32)bignum_digits(b)[0], vc;
unsigned32 vclow;
if (va < 0)
{ if (vb < 0) Dmultiply(vc, vclow, -va, -vb, 0);
else
{ Dmultiply(vc, vclow, -va, vb, 0);
if (vclow == 0) vc = -vc;
else
{ vclow = clear_top_bit(-(int32)vclow);
vc = ~vc;
}
}
}
else if (vb < 0)
{ Dmultiply(vc, vclow, va, -vb, 0);
if (vclow == 0) vc = -vc;
else
{ vclow = clear_top_bit(-(int32)vclow);
vc = ~vc;
}
}
else Dmultiply(vc, vclow, va, vb, 0);
return make_two_word_bignum(vc, vclow);
}
/*
* I take the absolute values of the two input values a and b,
* recording what the eventual sign for the product will need to be.
*/
if (((int32)bignum_digits(a)[lena-1]) < 0)
{ sign = -sign;
push(b);
/*
* Negating a negative bignum can sometimes mean that it will
* have to get longer (because of the twos complement assymmetry),
* but can never cause it to shrink, In particular it can never lead
* to demotion to a fixnum, so after this call to negateb it is still
* OK to assume that a is a bignum.
*/
a = negateb(a);
pop(b);
errexit();
lena = (bignum_length(a) - CELL)/4;
}
if (((int32)bignum_digits(b)[lenb-1]) < 0)
{ sign = -sign;
push(a);
/* see above comments about negateb */
b = negateb(b);
pop(a);
errexit();
lenb = (bignum_length(b) - CELL)/4;
}
if (lenb < lena) /* Commute so that b is at least as long as a */
{ c = a;
a = b;
b = c;
lenc = lena;
lena = lenb;
lenb = lenc;
}
push2(a, b);
/*
* For very big numbers I have two special actions called for here. First
* I must round up the size of the result vector to have a big enough power
* of two as a factor so that (recursive) splitting in two does not cause
* trouble later. Then I have to allocate some workspace, the size of that
* being related to the size of the numbers being handled.
*/
if (lena > KARATSUBA_CUTOFF)
{
int k = 0;
/*
* I pad lenc up to have a suitably large power of 2 as a factor so
* that splitting numbers in half works neatly for me.
*/
lenc = (lenb+1)/2; /* rounded up half-length of longer number */
while (lenc > KARATSUBA_CUTOFF)
{ lenc = (lenc + 1)/2;
k++;
}
while (k != 0)
{ lenc = 2*lenc;
k--;
}
lenc = 2*lenc;
c = getvector(TAG_NUMBERS, TYPE_BIGNUM, CELL+8*lenc);
errexitn(2);
/*
* The next line seems pretty shameless, but it may reduce the amount of
* garbage collection I do. When the workspace vector needed is short enough
* (at present up to 256 bytes) I use the character assembly buffer (boffo)
* as my workspace, relying on an expectation that bignumber multiplication
* can never be triggered from places within the reader where that buffer is
* in use for its normal purpose. I forge tag bits to make boffo look like
* a number here, but can never trigger garbage collection in a way that
* might inspect same, so that too is safe at present.
*/
if ((4*lenc+CELL) <= (intxx)length_of_header(vechdr(boffo)))
d = (Lisp_Object)((char *)boffo + TAG_NUMBERS - TAG_VECTOR);
else
{ push(c);
d = getvector(TAG_NUMBERS, TYPE_BIGNUM, CELL+4*lenc);
pop(c);
}
lenc = 2*lenc;
}
else
{
/*
* In cases where I will use classical long multiplication there is no
* need to waste space with extra padding or with the workspace vector d.
*/
lenc = lena + lenb;
c = getvector(TAG_NUMBERS, TYPE_BIGNUM, CELL+4*lenc);
d = c; /* set d to avoid dataflow anomaly */
}
pop2(b, a);
errexit();
{ unsigned32 *da = &bignum_digits(a)[0],
*db = &bignum_digits(b)[0],
*dc = &bignum_digits(c)[0],
*dd = &bignum_digits(d)[0];
long_times(dc, da, db, dd, lena, lenb, lenc);
}
/*
* Here the absolute value of the product has been computed properly.
* The result can easily have a zero top digit, which will need trimming
* off. If at least one of the input values was a number which had to
* be represented with a zero leading digit (e.g. 0x40000000) then the
* product can have two leading zero digits here. Similarly for negative
* results. Also padding (to allow splitting numbers into two) can have
* resulted in extra padding up at the top. lenc gives the size of vector
* that was allocated, lena+lenb is a much better guess of how much data
* is active in it.
*/
errexit();
{ int32 newlenc = lena + lenb;
/*
* I tidy up by putting a zero in any padding word above the top of the
* active data, and by inserting a header in space that gets trimmed off
* in such a way that the garbage collector will not get upset. This
* all just roughly trims the numbers - fine adjustment follows.
*/
#ifdef ADDRESS_64
if ((newlenc & 1) != 0)
#else
if ((newlenc & 1) == 0)
#endif
{ bignum_digits(c)[newlenc] = 0;
if (lenc != newlenc) /* i.e. I padded out somewhat */
{
*(Header *)&bignum_digits(c)[newlenc+1] =
make_bighdr(lenc-newlenc-1);
lenc = newlenc;
numhdr(c) = make_bighdr(lenc+CELL/4);
}
}
else if (lenc != newlenc) /* i.e. I padded out somewhat */
{
*(Header *)&bignum_digits(c)[newlenc] =
make_bighdr(lenc-newlenc);
lenc = newlenc;
numhdr(c) = make_bighdr(lenc+CELL/4);
}
}
/*
* Now I am safe against the garbage collector, and the number c has as
* its length just lena+lenb, even if it had been padded out earlier.
*/
if (sign < 0)
{ unsigned32 carry = 0xffffffffU;
for (i=0; i<lenc-1; i++)
{ carry = clear_top_bit(~bignum_digits(c)[i]) + top_bit(carry);
bignum_digits(c)[i] = clear_top_bit(carry);
}
carry = ~bignum_digits(c)[i] + top_bit(carry);
if (carry != -1)
{ bignum_digits(c)[i] = carry;
return c; /* no truncation needed */
}
carry = bignum_digits(c)[i-1];
if ((carry & 0x40000000) == 0)
{ bignum_digits(c)[i] = 0xffffffffU;
return c; /* no truncation becase of previous digit */
}
/*
* I need to argue that lenc was at least 2, so bignum_digits(c)[i-2]
* could at worst access the header word of the bignum - but it can never
* do that because if it were doing so then the bignum product would
* be about to have a value zero or thereabouts. One-word bignums are not
* allowed to have leading zero digits.
*/
if (carry == 0x7fffffff &&
(bignum_digits(c)[i-2] & 0x40000000) != 0) /* chop 2 */
{ bignum_digits(c)[i-2] |= ~0x7fffffff;
/*
* I common up the code to chop off two words from the number at label "chop2"
*/
goto chop2;
}
bignum_digits(c)[i-1] |= ~0x7fffffff;
/* Drop through to truncate by 1 and sometimes that is easy */
}
else
{ unsigned32 w = bignum_digits(c)[lenc-1];
if (w != 0) return c; /* no truncation */
w = bignum_digits(c)[lenc-2];
if ((w & 0x40000000) != 0) return c;
if (w == 0 &&
(bignum_digits(c)[lenc-3] & 0x40000000) == 0) goto chop2;
/* truncate one word */
}
/*
* here the data in the bignum is all correct (even in the most significant
* digit) but I need to shrink the number by one word. Because of all the
* doubleword alignment that is used here this can sometimes be done very
* easily, and other times it involves forging a short bit of dummy data
* to fill in a gap that gets left in the heap.
*/
numhdr(c) -= pack_hdrlength(1);
#ifdef ADDRESS_64
if ((lenc & 1) == 0)
#else
if ((lenc & 1) != 0)
#endif
bignum_digits(c)[lenc-1] = 0; /* tidy up */
else bignum_digits(c)[lenc-1] = make_bighdr(2);
return c;
chop2:
/*
* Trim two words from the number c
*/
numhdr(c) -= pack_hdrlength(2);
lenc -= 2;
bignum_digits(c)[lenc] = 0;
#ifdef ADDRESS_64
lenc = (lenc + 1) & ~1;
#else
lenc |= 1;
#endif
*(Header *)&bignum_digits(c)[lenc] = make_bighdr(2);
return c;
}
#ifdef COMMON
#define timesbr(a, b) timesir(a, b)
#define timesbc(a, b) timesic(a, b)
#endif
#define timesbf(a, b) timessf(a, b)
#ifdef COMMON
#define timesri(a, b) timesir(b, a)
#define timesrs(a, b) timessr(b, a)
#define timesrb(a, b) timesbr(b, a)
static Lisp_Object timesrr(Lisp_Object a, Lisp_Object b)
/*
* multiply a pair of rational numbers
*/
{
Lisp_Object nil = C_nil;
Lisp_Object w = nil;
push5(numerator(a), denominator(a),
numerator(b), denominator(b), 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;
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);
nil = C_nil;
if (exception_pending()) goto fail;
nb = quot2(nb, g);
nil = C_nil;
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);
fail:
popv(5);
return w;
#undef g
#undef db
#undef nb
#undef da
#undef na
}
#define timesrc(a, b) timesic(a, b)
#define timesrf(a, b) timessf(a, b)
#define timesci(a, b) timesic(b, a)
#define timescs(a, b) timessc(b, a)
#define timescb(a, b) timesbc(b, a)
#define timescr(a, b) timesrc(b, a)
static Lisp_Object timescc(Lisp_Object a, Lisp_Object b)
/*
* multiply a pair of complex values
*/
{
Lisp_Object nil = C_nil;
Lisp_Object w = nil;
push4(real_part(a), imag_part(a),
real_part(b), imag_part(b));
push2(nil, nil);
#define u stack[0]
#define v stack[-1]
#define ib stack[-2]
#define rb stack[-3]
#define ia stack[-4]
#define ra stack[-5]
u = times2(ra, rb);
nil = C_nil;
if (exception_pending()) goto fail;
v = times2(ia, ib);
nil = C_nil;
if (exception_pending()) goto fail;
v = negate(v);
nil = C_nil;
if (exception_pending()) goto fail;
u = plus2(u, v); /* real part of result */
nil = C_nil;
if (exception_pending()) goto fail;
v = times2(ra, ib);
nil = C_nil;
if (exception_pending()) goto fail;
ib = times2(rb, ia);
nil = C_nil;
if (exception_pending()) goto fail;
v = plus2(v, ib); /* imaginary part */
nil = C_nil;
if (exception_pending()) goto fail;
w = make_complex(u, v);
fail:
popv(6);
return w;
#undef u
#undef v
#undef ib
#undef rb
#undef ia
#undef ra
}
#define timescf(a, b) timesci(a, b)
#endif
#define timesfi(a, b) timesif(b, a)
#ifdef COMMON
#define timesfs(a, b) timessf(b, a)
#endif
#define timesfb(a, b) timesbf(b, a)
#ifdef COMMON
#define timesfr(a, b) timesrf(b, a)
#define timesfc(a, b) timescf(b, a)
#endif
static Lisp_Object timesff(Lisp_Object a, Lisp_Object b)
/*
* multiply boxed floats - see commentary on plusff()
*/
{
#ifdef COMMON
int32 ha = type_of_header(flthdr(a)), hb = type_of_header(flthdr(b));
#endif
double d = float_of_number(a) * float_of_number(b);
#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
}
/*
* ... and now for the dispatch code that copes with general
* multiplication.
*/
Lisp_Object times2(Lisp_Object a, Lisp_Object b)
{
switch ((int)a & TAG_BITS)
{
case TAG_FIXNUM:
switch ((int)b & TAG_BITS)
{
case TAG_FIXNUM:
return timesii(a, b);
#ifdef COMMON
case TAG_SFLOAT:
return timesis(a, b);
#endif
case TAG_NUMBERS:
{ int32 hb = type_of_header(numhdr(b));
switch (hb)
{
case TYPE_BIGNUM:
return timesib(a, b);
#ifdef COMMON
case TYPE_RATNUM:
return timesir(a, b);
case TYPE_COMPLEX_NUM:
return timesic(a, b);
#endif
default:
return aerror1("bad arg for times", b);
}
}
case TAG_BOXFLOAT:
return timesif(a, b);
default:
return aerror1("bad arg for times", b);
}
#ifdef COMMON
case TAG_SFLOAT:
switch (b & TAG_BITS)
{
case TAG_FIXNUM:
return timessi(a, b);
case TAG_SFLOAT:
{ Float_union aa, bb; /* timesss() coded in-line */
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 timessb(a, b);
case TYPE_RATNUM:
return timessr(a, b);
case TYPE_COMPLEX_NUM:
return timessc(a, b);
default:
return aerror1("bad arg for times", b);
}
}
case TAG_BOXFLOAT:
return timessf(a, b);
default:
return aerror1("bad arg for times", 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 timesbi(a, b);
#ifdef COMMON
case TAG_SFLOAT:
return timesbs(a, b);
#endif
case TAG_NUMBERS:
{ int32 hb = type_of_header(numhdr(b));
switch (hb)
{
case TYPE_BIGNUM:
return timesbb(a, b);
#ifdef COMMON
case TYPE_RATNUM:
return timesbr(a, b);
case TYPE_COMPLEX_NUM:
return timesbc(a, b);
#endif
default:
return aerror1("bad arg for times", b);
}
}
case TAG_BOXFLOAT:
return timesbf(a, b);
default:
return aerror1("bad arg for times", b);
}
#ifdef COMMON
case TYPE_RATNUM:
switch (b & TAG_BITS)
{
case TAG_FIXNUM:
return timesri(a, b);
case TAG_SFLOAT:
return timesrs(a, b);
case TAG_NUMBERS:
{ int32 hb = type_of_header(numhdr(b));
switch (hb)
{
case TYPE_BIGNUM:
return timesrb(a, b);
case TYPE_RATNUM:
return timesrr(a, b);
case TYPE_COMPLEX_NUM:
return timesrc(a, b);
default:
return aerror1("bad arg for times", b);
}
}
case TAG_BOXFLOAT:
return timesrf(a, b);
default:
return aerror1("bad arg for times", b);
}
case TYPE_COMPLEX_NUM:
switch (b & TAG_BITS)
{
case TAG_FIXNUM:
return timesci(a, b);
case TAG_SFLOAT:
return timescs(a, b);
case TAG_NUMBERS:
{ int32 hb = type_of_header(numhdr(b));
switch (hb)
{
case TYPE_BIGNUM:
return timescb(a, b);
case TYPE_RATNUM:
return timescr(a, b);
case TYPE_COMPLEX_NUM:
return timescc(a, b);
default:
return aerror1("bad arg for times", b);
}
}
case TAG_BOXFLOAT:
return timescf(a, b);
default:
return aerror1("bad arg for times", b);
}
#endif
default: return aerror1("bad arg for times", a);
}
}
case TAG_BOXFLOAT:
switch ((int)b & TAG_BITS)
{
case TAG_FIXNUM:
return timesfi(a, b);
#ifdef COMMON
case TAG_SFLOAT:
return timesfs(a, b);
#endif
case TAG_NUMBERS:
{ int32 hb = type_of_header(numhdr(b));
switch (hb)
{
case TYPE_BIGNUM:
return timesfb(a, b);
#ifdef COMMON
case TYPE_RATNUM:
return timesfr(a, b);
case TYPE_COMPLEX_NUM:
return timesfc(a, b);
#endif
default:
return aerror1("bad arg for times", b);
}
}
case TAG_BOXFLOAT:
return timesff(a, b);
default:
return aerror1("bad arg for times", b);
}
default:
return aerror1("bad arg for times", a);
}
}
/* end of arith02.c */