#include "Python.h"

#ifndef PY_NO_SHORT_FLOAT_REPR

#include "float.h"

#define MALLOC PyMem_Malloc

#define FREE PyMem_Free

#ifdef DOUBLE_IS_LITTLE_ENDIAN_IEEE754

# define IEEE_8087

#if defined(DOUBLE_IS_BIG_ENDIAN_IEEE754) || \

defined(DOUBLE_IS_ARM_MIXED_ENDIAN_IEEE754)

# define IEEE_MC68k

#if defined(IEEE_8087) + defined(IEEE_MC68k) != 1

#error "Exactly one of IEEE_8087 or IEEE_MC68k should be defined."

#if defined(WORDS_BIGENDIAN) && (defined(DOUBLE_IS_LITTLE_ENDIAN_IEEE754) || \

defined(DOUBLE_IS_ARM_MIXED_ENDIAN_IEEE754))

#error "doubles and ints have incompatible endianness"

#if !defined(WORDS_BIGENDIAN) && defined(DOUBLE_IS_BIG_ENDIAN_IEEE754)

#error "doubles and ints have incompatible endianness"

#if defined(HAVE_UINT32_T) && defined(HAVE_INT32_T)

typedef PY_UINT32_T ULong;

typedef PY_INT32_T Long;

#error "Failed to find an exact-width 32-bit integer type"

#if defined(HAVE_UINT64_T)

#define ULLong PY_UINT64_T

#undef ULLong

#undef DEBUG

#ifdef Py_DEBUG

#define DEBUG

#ifdef DEBUG

#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}

#ifndef PRIVATE_MEM

#define PRIVATE_MEM 2304

#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))

static double private_mem[PRIVATE_mem], *pmem_next = private_mem;

#ifdef __cplusplus

extern "C" {

typedef union { double d; ULong L[2]; } U;

#ifdef IEEE_8087

#define word0(x) (x)->L[1]

#define word1(x) (x)->L[0]

#define word0(x) (x)->L[0]

#define word1(x) (x)->L[1]

#define dval(x) (x)->d

#ifndef STRTOD_DIGLIM

#define STRTOD_DIGLIM 40

#ifndef MAX_ABS_EXP

#define MAX_ABS_EXP 19999U

#if defined(IEEE_8087)

#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \

((unsigned short *)a)[0] = (unsigned short)c, a++)

#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \

((unsigned short *)a)[1] = (unsigned short)c, a++)

#define Exp_shift 20

#define Exp_shift1 20

#define Exp_msk1 0x100000

#define Exp_msk11 0x100000

#define Exp_mask 0x7ff00000

#define P 53

#define Nbits 53

#define Bias 1023

#define Emax 1023

#define Emin (-1022)

#define Etiny (-1074) /* smallest denormal is 2**Etiny */

#define Exp_1 0x3ff00000

#define Exp_11 0x3ff00000

#define Ebits 11

#define Frac_mask 0xfffff

#define Frac_mask1 0xfffff

#define Ten_pmax 22

#define Bletch 0x10

#define Bndry_mask 0xfffff

#define Bndry_mask1 0xfffff

#define Sign_bit 0x80000000

#define Log2P 1

#define Tiny0 0

#define Tiny1 1

#define Quick_max 14

#define Int_max 14

#ifndef Flt_Rounds

#ifdef FLT_ROUNDS

#define Flt_Rounds FLT_ROUNDS

#define Flt_Rounds 1

#endif /*Flt_Rounds*/

#define Rounding Flt_Rounds

#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))

#define Big1 0xffffffff

typedef struct BCinfo BCinfo;

BCinfo {

int e0, nd, nd0, scale;

};

#define FFFFFFFF 0xffffffffUL

#define Kmax 7

Bigint {

};

typedef struct Bigint Bigint;

#ifndef Py_USING_MEMORY_DEBUGGER

static Bigint *freelist[Kmax+1];

static Bigint *

Balloc(int k)

{

int x;

Bigint *rv;

unsigned int len;

if (k <= Kmax && (rv = freelist[k]))

freelist[k] = rv->next;

else {

x = 1 << k;

len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)

/sizeof(double);

if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {

rv = (Bigint*)pmem_next;

pmem_next += len;

}

else {

rv = (Bigint*)MALLOC(len*sizeof(double));

if (rv == NULL)

return NULL;

}

rv->k = k;

rv->maxwds = x;

}

rv->sign = rv->wds = 0;

return rv;

}

static void

Bfree(Bigint *v)

{

if (v) {

if (v->k > Kmax)

FREE((void*)v);

else {

v->next = freelist[v->k];

freelist[v->k] = v;

}

}

}

static Bigint *

Balloc(int k)

{

int x;

Bigint *rv;

unsigned int len;

x = 1 << k;

len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)

/sizeof(double);

rv = (Bigint*)MALLOC(len*sizeof(double));

if (rv == NULL)

return NULL;

rv->k = k;

rv->maxwds = x;

rv->sign = rv->wds = 0;

return rv;

}

static void

Bfree(Bigint *v)

{

if (v) {

FREE((void*)v);

}

}

#endif /* Py_USING_MEMORY_DEBUGGER */

#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \

y->wds*sizeof(Long) + 2*sizeof(int))

static Bigint *

multadd(Bigint *b, int m, int a) /* multiply by m and add a */

{

int i, wds;

#ifdef ULLong

ULong *x;

ULLong carry, y;

ULong carry, *x, y;

ULong xi, z;

Bigint *b1;

wds = b->wds;

x = b->x;

i = 0;

carry = a;

do {

#ifdef ULLong

y = *x * (ULLong)m + carry;

carry = y >> 32;

*x++ = (ULong)(y & FFFFFFFF);

xi = *x;

y = (xi & 0xffff) * m + carry;

z = (xi >> 16) * m + (y >> 16);

carry = z >> 16;

*x++ = (z << 16) + (y & 0xffff);

}

while(++i < wds);

if (carry) {

if (wds >= b->maxwds) {

b1 = Balloc(b->k+1);

if (b1 == NULL){

Bfree(b);

return NULL;

}

Bcopy(b1, b);

Bfree(b);

b = b1;

}

b->x[wds++] = (ULong)carry;

b->wds = wds;

}

return b;

}

static Bigint *

s2b(const char *s, int nd0, int nd, ULong y9)

{

Bigint *b;

int i, k;

Long x, y;

x = (nd + 8) / 9;

for(k = 0, y = 1; x > y; y <<= 1, k++) ;

b = Balloc(k);

if (b == NULL)

return NULL;

b->x[0] = y9;

b->wds = 1;

if (nd <= 9)

return b;

s += 9;

for (i = 9; i < nd0; i++) {

b = multadd(b, 10, *s++ - '0');

if (b == NULL)

return NULL;

}

s++;

for(; i < nd; i++) {

b = multadd(b, 10, *s++ - '0');

if (b == NULL)

return NULL;

}

return b;

}

static int

hi0bits(ULong x)

{

int k = 0;

if (!(x & 0xffff0000)) {

k = 16;

x <<= 16;

}

if (!(x & 0xff000000)) {

k += 8;

x <<= 8;

}

if (!(x & 0xf0000000)) {

k += 4;

x <<= 4;

}

if (!(x & 0xc0000000)) {

k += 2;

x <<= 2;

}

if (!(x & 0x80000000)) {

k++;

if (!(x & 0x40000000))

return 32;

}

return k;

}

static int

lo0bits(ULong *y)

{

int k;

ULong x = *y;

if (x & 7) {

if (x & 1)

return 0;

if (x & 2) {

*y = x >> 1;

return 1;

}

*y = x >> 2;

return 2;

}

k = 0;

if (!(x & 0xffff)) {

k = 16;

x >>= 16;

}

if (!(x & 0xff)) {

k += 8;

x >>= 8;

}

if (!(x & 0xf)) {

k += 4;

x >>= 4;

}

if (!(x & 0x3)) {

k += 2;

x >>= 2;

}

if (!(x & 1)) {

k++;

x >>= 1;

if (!x)

return 32;

}

*y = x;

return k;

}

static Bigint *

i2b(int i)

{

Bigint *b;

b = Balloc(1);

if (b == NULL)

return NULL;

b->x[0] = i;

b->wds = 1;

return b;

}

static Bigint *

mult(Bigint *a, Bigint *b)

{

Bigint *c;

int k, wa, wb, wc;

ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;

ULong y;

#ifdef ULLong

ULLong carry, z;

ULong carry, z;

ULong z2;

if ((!a->x[0] && a->wds == 1) || (!b->x[0] && b->wds == 1)) {

c = Balloc(0);

if (c == NULL)

return NULL;

c->wds = 1;

c->x[0] = 0;

return c;

}

if (a->wds < b->wds) {

c = a;

a = b;

b = c;

}

k = a->k;

wa = a->wds;

wb = b->wds;

wc = wa + wb;

if (wc > a->maxwds)

k++;

c = Balloc(k);

if (c == NULL)

return NULL;

for(x = c->x, xa = x + wc; x < xa; x++)

*x = 0;

xa = a->x;

xae = xa + wa;

xb = b->x;

xbe = xb + wb;

xc0 = c->x;

#ifdef ULLong

for(; xb < xbe; xc0++) {

if ((y = *xb++)) {

x = xa;

xc = xc0;

carry = 0;

do {

z = *x++ * (ULLong)y + *xc + carry;

carry = z >> 32;

*xc++ = (ULong)(z & FFFFFFFF);

}

while(x < xae);

*xc = (ULong)carry;

}

}

for(; xb < xbe; xb++, xc0++) {

if (y = *xb & 0xffff) {

x = xa;

xc = xc0;

carry = 0;

do {

z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;

carry = z >> 16;

z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;

carry = z2 >> 16;

Storeinc(xc, z2, z);

}

while(x < xae);

*xc = carry;

}

if (y = *xb >> 16) {

x = xa;

xc = xc0;

carry = 0;

z2 = *xc;

do {

z = (*x & 0xffff) * y + (*xc >> 16) + carry;

carry = z >> 16;

Storeinc(xc, z, z2);

z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;

carry = z2 >> 16;

}

while(x < xae);

*xc = z2;

}

}

for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;

c->wds = wc;

return c;

}

#ifndef Py_USING_MEMORY_DEBUGGER

static Bigint *p5s;

static Bigint *

pow5mult(Bigint *b, int k)

{

Bigint *b1, *p5, *p51;

int i;

static int p05[3] = { 5, 25, 125 };

if ((i = k & 3)) {

b = multadd(b, p05[i-1], 0);

if (b == NULL)

return NULL;

}

if (!(k >>= 2))

return b;

p5 = p5s;

if (!p5) {

/* first time */

p5 = i2b(625);

if (p5 == NULL) {

Bfree(b);

return NULL;

}

p5s = p5;

p5->next = 0;

}

for(;;) {

if (k & 1) {

b1 = mult(b, p5);

Bfree(b);

b = b1;

if (b == NULL)

return NULL;

}

if (!(k >>= 1))

break;

p51 = p5->next;

if (!p51) {

p51 = mult(p5,p5);

if (p51 == NULL) {

Bfree(b);

return NULL;

}

p51->next = 0;

p5->next = p51;

}

p5 = p51;

}

return b;

}

static Bigint *

pow5mult(Bigint *b, int k)

{

Bigint *b1, *p5, *p51;

int i;

static int p05[3] = { 5, 25, 125 };

if ((i = k & 3)) {

b = multadd(b, p05[i-1], 0);

if (b == NULL)

return NULL;

}

if (!(k >>= 2))

return b;

p5 = i2b(625);

if (p5 == NULL) {

Bfree(b);

return NULL;

}

for(;;) {

if (k & 1) {

b1 = mult(b, p5);

Bfree(b);

b = b1;

if (b == NULL) {

Bfree(p5);

return NULL;

}

}

if (!(k >>= 1))

break;

p51 = mult(p5, p5);

Bfree(p5);

p5 = p51;

if (p5 == NULL) {

Bfree(b);

return NULL;

}

}

Bfree(p5);

return b;

}

#endif /* Py_USING_MEMORY_DEBUGGER */

static Bigint *

lshift(Bigint *b, int k)

{

int i, k1, n, n1;

Bigint *b1;

ULong *x, *x1, *xe, z;

if (!k || (!b->x[0] && b->wds == 1))

return b;

n = k >> 5;

k1 = b->k;

n1 = n + b->wds + 1;

for(i = b->maxwds; n1 > i; i <<= 1)

k1++;

b1 = Balloc(k1);

if (b1 == NULL) {

Bfree(b);

return NULL;

}

x1 = b1->x;

for(i = 0; i < n; i++)

*x1++ = 0;

x = b->x;

xe = x + b->wds;

if (k &= 0x1f) {

k1 = 32 - k;

z = 0;

do {

*x1++ = *x << k | z;

z = *x++ >> k1;

}

while(x < xe);

if ((*x1 = z))

++n1;

}

else do

*x1++ = *x++;

while(x < xe);

b1->wds = n1 - 1;

Bfree(b);

return b1;

}

static int

cmp(Bigint *a, Bigint *b)

{

ULong *xa, *xa0, *xb, *xb0;

int i, j;

i = a->wds;

j = b->wds;

#ifdef DEBUG

if (i > 1 && !a->x[i-1])

Bug("cmp called with a->x[a->wds-1] == 0");

if (j > 1 && !b->x[j-1])

Bug("cmp called with b->x[b->wds-1] == 0");

if (i -= j)

return i;

xa0 = a->x;

xa = xa0 + j;

xb0 = b->x;

xb = xb0 + j;

for(;;) {

if (*--xa != *--xb)

return *xa < *xb ? -1 : 1;

if (xa <= xa0)

break;

}

return 0;

}

static Bigint *

diff(Bigint *a, Bigint *b)

{

Bigint *c;

int i, wa, wb;

ULong *xa, *xae, *xb, *xbe, *xc;

#ifdef ULLong

ULLong borrow, y;

ULong borrow, y;

ULong z;

i = cmp(a,b);

if (!i) {

c = Balloc(0);

if (c == NULL)

return NULL;

c->wds = 1;

c->x[0] = 0;

return c;

}

if (i < 0) {

c = a;

a = b;

b = c;

i = 1;

}

else

i = 0;

c = Balloc(a->k);

if (c == NULL)

return NULL;

c->sign = i;

wa = a->wds;

xa = a->x;

xae = xa + wa;

wb = b->wds;

xb = b->x;

xbe = xb + wb;

xc = c->x;

borrow = 0;

#ifdef ULLong

do {

y = (ULLong)*xa++ - *xb++ - borrow;

borrow = y >> 32 & (ULong)1;

*xc++ = (ULong)(y & FFFFFFFF);

}

while(xb < xbe);

while(xa < xae) {

y = *xa++ - borrow;

borrow = y >> 32 & (ULong)1;

*xc++ = (ULong)(y & FFFFFFFF);

}

do {

y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;

borrow = (y & 0x10000) >> 16;

z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;

borrow = (z & 0x10000) >> 16;

Storeinc(xc, z, y);

}

while(xb < xbe);

while(xa < xae) {

y = (*xa & 0xffff) - borrow;

borrow = (y & 0x10000) >> 16;

z = (*xa++ >> 16) - borrow;

borrow = (z & 0x10000) >> 16;

Storeinc(xc, z, y);

}

while(!*--xc)

wa--;

c->wds = wa;

return c;

}

static double

ulp(U *x)

{

Long L;

U u;

L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;

word0(&u) = L;

word1(&u) = 0;

return dval(&u);

}

static double

b2d(Bigint *a, int *e)