mirror of
https://gitlab.com/embeddable-common-lisp/ecl.git
synced 2026-01-02 07:30:55 -08:00
494 lines
9.8 KiB
C
494 lines
9.8 KiB
C
/* -*- mode: c; c-basic-offset: 8 -*- */
|
|
/*
|
|
earith.c -- Support for bignum arithmetic.
|
|
*/
|
|
/*
|
|
Copyright (c) 1984, Taiichi Yuasa and Masami Hagiya.
|
|
Copyright (c) 1990, Giuseppe Attardi.
|
|
Copyright (c) 2001, Juan Jose Garcia Ripoll.
|
|
|
|
ECL is free software; you can redistribute it and/or
|
|
modify it under the terms of the GNU Library General Public
|
|
License as published by the Free Software Foundation; either
|
|
version 2 of the License, or (at your option) any later version.
|
|
|
|
See file '../Copyright' for full details.
|
|
*/
|
|
|
|
/*
|
|
|
|
EXTENDED_MUL and EXTENDED_DIV perform 32 bit multiplication and
|
|
division, respectively.
|
|
|
|
EXTENDED_MUL(D,Q,R,HP,LP)
|
|
calculates D*Q+R and saves the result into the locations HP and LP.
|
|
D, Q, and R are 32 bit non-negative integers and HP and LP are
|
|
word addresses. The word at LP will contain the lower 31 (not 32)
|
|
bits of the result and its most significant bit is set 0. The word
|
|
at HP will contain the rest of the result and its MSB is also set 0.
|
|
|
|
EXTENDED_DIV(D,H,L,QP,RP)
|
|
divides [H:L] by D and saves the quotient and the remainder into
|
|
the locations QP and RP, respectively. D, H, L are 32 bit non-negative
|
|
integers and QP and RP are word addresses. Here, [H:L] means the
|
|
64 bit integer (imaginary) represented by H and L as follows.
|
|
|
|
63 62 31 30 0
|
|
|0|0|<lower 31 bits of H>|<lower 31 bits of L>|
|
|
|
|
Although [H:L] is 64 bits, you can assume that the quotient is always
|
|
represented as 32 bit non-negative integer.
|
|
*/
|
|
|
|
#include <ecl/ecl.h>
|
|
|
|
#ifdef CONVEX
|
|
|
|
static void
|
|
extended_mul(int d, int q, int r, int *hp, int *lp)
|
|
{
|
|
long long int ld, lq, lr, z;
|
|
int zh, zl;
|
|
|
|
ld = d;
|
|
lq = q;
|
|
lr = r;
|
|
z = ld*lq+lr;
|
|
zl = (z & 0x000000007fffffffLL);
|
|
zh = (z >> 31LL);
|
|
*hp = zh;
|
|
*lp = zl;
|
|
}
|
|
|
|
static void
|
|
extended_div(int d, int h, int l, int *qp, int *rp)
|
|
{
|
|
long long int lh, ld, ll;
|
|
|
|
ld = d;
|
|
lh = h;
|
|
ll = l;
|
|
lh = (lh << 31LL);
|
|
lh = (lh | ll);
|
|
*qp = (lh/ld);
|
|
*rp = (lh%ld);
|
|
}
|
|
#endif /* CONVEX */
|
|
|
|
#ifdef i386
|
|
|
|
static void
|
|
extended_mul(int d, int q, int r, int *hp, int *lp)
|
|
{ asm("pushl %ecx");
|
|
asm("movl 8(%ebp),%eax");
|
|
asm("mull 12(%ebp)");
|
|
asm("addl 16(%ebp),%eax");
|
|
asm("adcl $0,%edx");
|
|
asm("shll $1,%edx");
|
|
asm("btrl $31,%eax");
|
|
asm("adcl $0,%edx");
|
|
asm("movl 20(%ebp),%ecx");
|
|
asm("movl %edx, (%ecx)");
|
|
asm("movl 24(%ebp), %ecx");
|
|
asm("movl %eax, (%ecx)");
|
|
asm("popl %ecx");
|
|
}
|
|
|
|
static void
|
|
extended_div(int d, int h, int l, int *qp, int *rp)
|
|
{
|
|
asm("pushl %ebx");
|
|
asm("movl 12(%ebp),%edx");
|
|
asm("movl 16(%ebp),%eax");
|
|
asm("btl $0,%edx");
|
|
asm("jae 1f");
|
|
asm("btsl $31,%eax");
|
|
asm("1: shrl $1,%edx");
|
|
asm("idivl 8(%ebp)");
|
|
asm("movl 20(%ebp),%ebx");
|
|
asm("movl %eax,(%ebx)");
|
|
asm("movl 24(%ebp),%ebx");
|
|
asm("movl %edx,(%ebx)");
|
|
asm("popl %ebx");
|
|
}
|
|
#endif /* i386 */
|
|
|
|
#ifdef IBMRT
|
|
|
|
static void
|
|
extended_mul(int d, int q, int r, int *hp, int *lp)
|
|
{
|
|
/* d=L750+20, q=L750+24, etc. */
|
|
asm(" get r0,L750+20(r13)"); /* get an argument */
|
|
asm(" mts r10,r0"); /* put in MQ */
|
|
asm(" get r2,L750+24(r13)"); /* get the other argument */
|
|
asm(" s r0,r0"); /* zero partial product. set carry to 1. */
|
|
asm(" m r0,r2
|
|
m r0,r2
|
|
m r0,r2
|
|
m r0,r2
|
|
m r0,r2
|
|
m r0,r2
|
|
m r0,r2
|
|
m r0,r2
|
|
m r0,r2
|
|
m r0,r2
|
|
m r0,r2
|
|
m r0,r2
|
|
m r0,r2
|
|
m r0,r2
|
|
m r0,r2
|
|
m r0,r2");
|
|
/* Now (r0)//mq has the 64 bit product; overflow is ignored. */
|
|
asm(" mfs r10,r2"); /* set r2 = low order word of result
|
|
* so product is in (r0)//(r2).
|
|
*/
|
|
/*
|
|
* Force product into two single precision words.
|
|
*/
|
|
|
|
asm(" get r3,$1f
|
|
sli r0,1
|
|
ar2,r2
|
|
bnc0r r3"); /* branch if carry = 0 */
|
|
asm(" oil r0,r0,1
|
|
1:
|
|
sri r2,1");
|
|
/* Now add in the third argument. */
|
|
asm(" get r4,$2f
|
|
get r3,L750+28(r13)
|
|
a r2,r3
|
|
|
|
bnmr r4"); /* branch if not minus */
|
|
asm(" clrbu r2,0
|
|
lis r3,1
|
|
a r0,r3
|
|
2:
|
|
|
|
get r3,L750+32(r13)
|
|
put r0,0(r3)
|
|
get r3,L750+36(r13)
|
|
put r2,0(r3)
|
|
");
|
|
}
|
|
|
|
static void
|
|
extended_div(int d, int h, int l, int *qp, int *rp)
|
|
{
|
|
/* d=L754+20, h=L754+24, etc. */
|
|
/* Move arguments into registers. */
|
|
asm(" get r0,L754+28(r13)"); /* Low order word of dividend. */
|
|
asm(" get r2,L754+24(r13)"); /* High order word of dividend. */
|
|
asm(" mttbil r2,15
|
|
mftbiu r0,0
|
|
sri r2,1
|
|
mts r10,r0
|
|
get r3,L754+20(r13)") /* Divisor. */
|
|
/* Perform 32 steps of division. */
|
|
asm(" d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3
|
|
d r2,r3");
|
|
/* Now MQ has the quotient, R2 the remainder, and R3 is
|
|
* the unchanged divisor. */
|
|
asm(" mttbiu r2,0"); /* Do add-back if necessary. */
|
|
asm(" jntb 1f
|
|
a r2,r3
|
|
1:
|
|
mfs r10,r0
|
|
c r2,r3"); /* Remainder : divisor. */
|
|
asm(" jne 2f
|
|
inc r0,1
|
|
x r2,r2
|
|
2:");
|
|
/* Now r0 has the quotient and r2 has the remainder. */
|
|
asm(" get r3,L754+32(r13)"); /* Quotient address. */
|
|
asm(" put r0,0(r3)");
|
|
asm(" get r3,L754+36(r13)"); /* Remainder address. */
|
|
asm(" put r2,0(r3)");
|
|
}
|
|
|
|
#endif /* IBMRT */
|
|
|
|
#if defined(NEWS) || defined(MAC)
|
|
|
|
static void
|
|
extended_mul(int d, int q, int r, int *hp, int *lp)
|
|
{
|
|
asm(" move.l d2,-(sp)
|
|
clr.l d2
|
|
move.l (8,fp),d0
|
|
mulu.l (12,fp),d1:d0
|
|
add.l (16,fp),d0
|
|
addx.l d2,d1
|
|
lsl.l #1,d0
|
|
roxl.l #1,d1
|
|
lsr.l #1,d0
|
|
move.l (20,fp),a0
|
|
move.l d1,(a0)
|
|
move.l (24,a6),a0
|
|
move.l d0,(a0)");
|
|
}
|
|
|
|
static void
|
|
extended_div(int d, int h, int l, int *qp, int *rp)
|
|
{
|
|
asm("movem.l (12,fp),#0x303
|
|
lsl.l #1,d1
|
|
lsr.l #1,d0
|
|
roxr.l #1,d1
|
|
divu.l (8,fp),d0:d1
|
|
move.l d0,(a1)
|
|
move.l d1,(a0)
|
|
");
|
|
}
|
|
|
|
#endif /* NEWS || MAC */
|
|
|
|
#ifdef __mips
|
|
|
|
/* earith.s for MIPS R2000 processor
|
|
by Doug Katzman
|
|
version 2.1.d dated 7/13/89 15:31 EDT */
|
|
|
|
/* Register names:
|
|
#define v0 $2 return value
|
|
#define v1 $3
|
|
#define a0 $4 argument registers
|
|
#define a1 $5
|
|
#define a2 $6
|
|
#define a3 $7
|
|
#define t7 $15
|
|
*/
|
|
|
|
static void
|
|
extended_mul(unsigned int d, unsigned int q, unsigned int r, unsigned int *hp,
|
|
unsigned int *lp)
|
|
{
|
|
asm("mult $4, $5"); /* [hi:lo] = d * q */
|
|
asm("mfhi $5"); /* a1 = hi */
|
|
asm("sll $5, 1");
|
|
asm("mflo $4");
|
|
asm("srl $15, $4, 31");
|
|
asm("and $4, 0x7fffffff");
|
|
asm("or $5, $15");
|
|
asm("addu $4, $6"); /* [a1:a0] += r */
|
|
asm("srl $15, $4, 31");
|
|
asm("and $4, 0x7fffffff");
|
|
asm("addu $5, $15");
|
|
asm("sw $5, 0($7)"); /* *hp = a1 */
|
|
#ifdef __GNUC__
|
|
asm("lw $7, %0" :: "g" (lp));
|
|
#else
|
|
asm("lw $7, 16($sp)"); /* fetch fifth actual argument from stack */
|
|
#endif
|
|
asm("sw $4, 0($7)"); /* *lp = a0 */
|
|
}
|
|
|
|
static void
|
|
extended_div(unsigned int d, unsigned int h, unsigned int l, unsigned int *qp,
|
|
unsigned int *rp)
|
|
{
|
|
|
|
asm("sll $6, 1");
|
|
asm("li $2, 31"); /* v0 holds number of shifts */
|
|
asm("loop:
|
|
srl $15, $6, 31");
|
|
asm("sll $5, 1");
|
|
asm("or $5, $15");
|
|
asm("sll $6, 1");
|
|
asm("subu $15, $5, $4"); /* t = h - d */
|
|
asm("bltz $15, underflow");
|
|
asm("move $5, $15");
|
|
asm("or $6, 1");
|
|
asm("underflow:
|
|
subu $2, 1");
|
|
asm("bnez $2, loop");
|
|
asm("sw $6, 0($7)"); /* *qp = l */
|
|
#ifdef __GNUC__
|
|
asm("lw $7, %0" :: "g" (rp));
|
|
#else
|
|
asm("lw $7, 16($sp)"); /* fetch fifth actual argument from stack */
|
|
#endif
|
|
asm("sw $5, 0($7)"); /* *rp = h */
|
|
}
|
|
#endif /* __mips */
|
|
|
|
#if defined(sun3) || (defined __NeXT)
|
|
|
|
static void
|
|
extended_mul(int d, int q, int r, int *hp, int *lp)
|
|
{
|
|
asm("
|
|
movl d2,a7@-
|
|
clrl d2
|
|
movl a6@(8),d0
|
|
mulul a6@(12),d1:d0
|
|
addl a6@(16),d0
|
|
addxl d2,d1
|
|
lsll #1,d0
|
|
roxll #1,d1
|
|
lsrl #1,d0
|
|
movl a6@(20),a0
|
|
movl d1,a0@
|
|
movl a6@(24),a0
|
|
movl d0,a0@
|
|
movl a7@+,d2
|
|
");
|
|
}
|
|
|
|
static void
|
|
extended_div(int d, int h, int l, int *qp, int *rp)
|
|
{
|
|
asm("moveml a6@(12),#0x303
|
|
lsll #1,d1
|
|
lsrl #1,d0
|
|
roxrl #1,d1
|
|
divul a6@(8),d0:d1
|
|
movl d0,a1@
|
|
movl d1,a0@
|
|
");
|
|
}
|
|
|
|
#endif /* sun3 */
|
|
|
|
/* Possible assembler version:
|
|
#ifdef sparc
|
|
_extended_mul:
|
|
!#PROLOGUE# 0
|
|
!#PROLOGUE# 1
|
|
save %sp,-96,%sp
|
|
mov %i0,%o0
|
|
call .umul,2
|
|
mov %i1,%o1
|
|
addcc %o0,%i2,%i0
|
|
addx %o1,0,%o1
|
|
sll %o1,1,%o1
|
|
tst %i0
|
|
bge L77003
|
|
sethi %hi(0x7fffffff),%o3
|
|
or %o3,%lo(0x7fffffff),%o3 ! [internal]
|
|
and %i0,%o3,%i0
|
|
inc %o1
|
|
L77003:
|
|
st %i0,[%i4]
|
|
st %o1,[%i3]
|
|
ret
|
|
restore %g0,0,%o0
|
|
|
|
#endif sparc
|
|
*/
|
|
|
|
#if defined(sparc) || defined(APOLLO) || defined(hpux) || defined(UNIGRAPH)n
|
|
|
|
/* for the time being use the C version:*/
|
|
|
|
static void
|
|
extended_mul(unsigned int d, unsigned int q, unsigned int r, unsigned int *hp,
|
|
unsigned int *lp)
|
|
{
|
|
register unsigned short dlo = d & 0xffff,
|
|
dhi = d >> 16,
|
|
qlo = q & 0xffff,
|
|
qhi = q >> 16;
|
|
unsigned int d0 = dhi * qlo + dlo * qhi,
|
|
d1 = dhi * qhi,
|
|
d2 = dlo * qlo;
|
|
|
|
d1 = (d1 << 1) + (d0 >> 15); /* add 17 MSB of d0 */
|
|
d1 += d2 >> 31; /* add MSB of d2 */
|
|
d2 &= 0x7fffffff; /* clear MSB of d2 */
|
|
d2 += (d0 & 0x7fff) << 16; /* add 15 LSB of d0: no overflow occurs */
|
|
d1 += d2 >> 31; /* add MSB of d2 */
|
|
d2 &= 0x7fffffff; /* clear MSB of d2 */
|
|
d2 += r;
|
|
d1 += d2 >> 31; /* add MSB of d2 */
|
|
d2 &= 0x7fffffff; /* clear MSB of d2 */
|
|
|
|
*hp = d1;
|
|
*lp = d2;
|
|
}
|
|
|
|
static void
|
|
extended_div(unsigned int d, unsigned int h, unsigned int l, unsigned int *qp,
|
|
unsigned int *rp)
|
|
{
|
|
int i;
|
|
int borrow;
|
|
|
|
l = (l << 1) | 1;
|
|
|
|
for (i = 31; i >= 0;) {
|
|
|
|
if (h >= d) {
|
|
h -= d;
|
|
borrow = 0;
|
|
}
|
|
else
|
|
borrow = 1;
|
|
|
|
if (i--)
|
|
h = (h << 1) | ((unsigned)l >> 31);
|
|
|
|
l = (l << 1) | borrow;
|
|
|
|
}
|
|
|
|
l = - l - 1;
|
|
|
|
*qp = l;
|
|
*rp = h;
|
|
}
|
|
|
|
#endif /* sparc */
|
|
|
|
#ifdef vax
|
|
|
|
static void
|
|
extended_mul(int d, int q, int r, int *hp, int *lp)
|
|
{
|
|
asm(" emul 4(ap),8(ap),12(ap),r0");
|
|
asm(" ashq $1,r0,r0");
|
|
asm(" rotl $-1,r0,r0");
|
|
asm(" movl r0,*20(ap)");
|
|
asm(" movl r1,*16(ap)");
|
|
}
|
|
|
|
static void
|
|
extended_div(int d, int h, int l, int *qp, int *rp)
|
|
{
|
|
asm(" clrl r0");
|
|
asm(" movl 8(ap),r1");
|
|
asm(" ashq $-1,r0,r0");
|
|
asm(" addl2 12(ap),r0");
|
|
asm(" ediv 4(ap),r0,*16(ap),*20(ap)");
|
|
}
|
|
#endif /* vax */
|