/*-----------------------------------------------------------------*-C-*---
 * File:    handc/runtime/basemath.c
 *
 *          Copyright (C)1997 Donovan Kolbly <d.kolbly@rscheme.org>
 *          as part of the RScheme project, licensed for free use.
 *          See <http://www.rscheme.org/> for the latest information.
 *
 * File version:     1.7
 * File mod date:    1997.11.29 23:10:50
 * System build:     v0.7.2, 97.12.21
 *
 * Purpose:          Arithmetic for basic (built-in) number system
 *------------------------------------------------------------------------*/

#include <rscheme/scheme.h>
#include <rscheme/smemory.h>

#if !FULL_NUMERIC_TOWER

void init_math( void )
{
}

#include "numsimpl.ci"

#else
#include <gmp.h>

typedef UINT_32  UWtype;   /* define types that longlong.h needs */


typedef UINT_16  UHWtype;
typedef INT_32   SItype;
typedef UINT_32  USItype;

#define W_TYPE_SIZE   (32)

#include <longlong.h>   /* from the GMP distribution */
/***********************************************************************
 Note:  There is a bug in longlong.h from the 2.0.2 version of GMP.
        Consider applying the following patch:

*** longlong.h  Sat Oct 25 22:47:20 1997
--- /tmp/gmp-2.0.2.orig/longlong.h      Fri May 24 01:34:24 1996
***************
*** 1295 ****
! #if !defined (smul_ppmm)
--- 1295 ----
! #if !defined (umul_ppmm)

***********************************************************************/

#define BIGNUM_P(x)    OBJ_ISA_PTR_OF_CLASS(x,bignum_class)
#define RATIONAL_P(x)  OBJ_ISA_PTR_OF_CLASS(x,mp_rational_class)
#define COMPLEX_P(x)   OBJ_ISA_PTR_OF_CLASS(x,rect_complex_class)

static void *mp_alloc( size_t n )
{
  return PTR_TO_DATAPTR(alloc( n, mp_data_class ));
}

static void *mp_realloc( void *old_data, size_t old_size, size_t new_size )
{
  void *new_data = PTR_TO_DATAPTR(alloc( new_size, mp_data_class ));
  memcpy( new_data, old_data, old_size );
  return new_data;
}

static void mp_free( void *ptr, size_t n )
{
  /* do nothing -- data will get GC'd */
}

void init_math( void )
{
  mp_set_memory_functions( mp_alloc, mp_realloc, mp_free );
}

#define DATAPTR_TO_PTR(ptr) OBJ(POINTER_TAG+(UINT_32)ptr)

static obj mpz_to_bignum( mpz_t n )
{
  return make3( bignum_class, 
		int2fx( n[0]._mp_alloc ),
		int2fx( n[0]._mp_size ),
		DATAPTR_TO_PTR( n[0]._mp_d ) );
}


static obj fx_to_bignum( obj x )
{
  mpz_t n;
  mpz_init_set_si( n, fx2int(x) );
  return mpz_to_bignum( n );
}

#define longint_to_float(l) int_64_to_float(extract_int_64(l))
#define longint_to_rational(l) bignum_to_rational(longint_to_bignum(l))

static obj longint_to_bignum( obj lng )
{
  return fx_to_bignum(ZERO);
}

static obj bignum_to_rational( obj b )
{
  return make2( mp_rational_class, b, fx_to_bignum( int2fx(1) ) );
}

static obj fx_to_rational( obj fx )
{
  return bignum_to_rational( fx_to_bignum(fx) );
}

#define STUBBED_OUT(r)  scheme_error( "basemath.c:~d: function stubbed out", 1, int2fx(__LINE__) ); return r


/************************ INT_32 OPERATIONS ************************/

#define HIGHBIT (1<<29)

static _rs_inline obj int_plus( INT_32 a, INT_32 b )
{
  INT_32 c = a + b;
  if ((~(a ^ b) & (c ^ a)) & HIGHBIT)
    {
      INT_64 a2 = int_32_to_int_64(a);
      INT_64 b2 = int_32_to_int_64(b);
      return int_64_compact( int_64_add( a2, b2 ) );
    }
  return int2fx( c );
}

static _rs_inline obj int_minus( INT_32 a, INT_32 b )
{
  INT_32 c = a - b;
  if ((~(a ^ b) & (c ^ a)) & HIGHBIT)
    {
      INT_64 a2 = int_32_to_int_64(a);
      INT_64 b2 = int_32_to_int_64(b);
      return int_64_compact( int_64_sub( a2, b2 ) );
    }
  return int2fx( a-b );
}

static _rs_inline obj int_mul( INT_32 a, INT_32 b )
{
  INT_32 p_hi, p_lo;
  INT_64 a2, b2;

  smul_ppmm( p_hi, p_lo, a, b );
  if (p_hi == 0)
    {
      if (p_lo < HIGHBIT)
	{
	  return int2fx( p_lo );
	}
    }
  else if (p_hi == -1)
    {
      if (p_lo >= -HIGHBIT)
	{
	  return int2fx( p_lo );
	}
    }
  a2 = int_32_to_int_64(a);
  b2 = int_32_to_int_64(b);
  return int_64_compact( int_64_mul( a2, b2 ) );
}

static _rs_inline obj int_div( INT_32 a, INT_32 b )
{
  return make_float( (IEEE_64)a/(IEEE_64)b );
}

static _rs_inline int int_cmp( INT_32 a, INT_32 b )
{
  if (a < b)
    return -1;
  else if (a > b)
    return 1;
  else
    return 0;
}

/************************ INT_64 OPERATIONS ************************/

static _rs_inline obj long_plus( INT_64 a, INT_64 b )
{
  return int_64_compact( int_64_add(a,b) );
}

static _rs_inline obj long_minus( INT_64 a, INT_64 b )
{
  return int_64_compact( int_64_sub(a,b) );
}

static _rs_inline obj long_mul( INT_64 a, INT_64 b )
{
  return int_64_compact( int_64_mul(a,b) );
}

static _rs_inline obj long_div( INT_64 a, INT_64 b )
{
  return make_float( int_64_to_float(a) / int_64_to_float(b) );
}

#define long_cmp(a,b) int_64_cmp(a,b)

/************************ IEEE_64 OPERATIONS ************************/

static _rs_inline obj fl_plus( IEEE_64 a, IEEE_64 b )
{
  return make_float( a + b );
}

static _rs_inline obj fl_minus( IEEE_64 a, IEEE_64 b )
{
  return make_float( a - b );
}

static _rs_inline obj fl_mul( IEEE_64 a, IEEE_64 b )
{
  return make_float( a * b );
}

static _rs_inline obj fl_div( IEEE_64 a, IEEE_64 b )
{
  return make_float( a / b );
}

static _rs_inline int fl_cmp( IEEE_64 a, IEEE_64 b )
{
  if (a < b)
    return -1;
  else if (a > b)
    return 1;
  else
    return 0;
}

/************************ BIGNUM OPERATIONS ************************/

static _rs_inline obj bignum_plus( obj a, obj b )
{
  STUBBED_OUT(ZERO);
}

static _rs_inline obj bignum_minus( obj a, obj b )
{
  STUBBED_OUT(ZERO);
}

static _rs_inline obj bignum_mul( obj a, obj b )
{
  STUBBED_OUT(ZERO);
}

static _rs_inline obj bignum_div( obj a, obj b )
{
  STUBBED_OUT(ZERO);
}

static int bignum_cmp( obj a, obj b )
{
  STUBBED_OUT(0);
}

/************************ RATIONAL OPERATIONS ************************/

static _rs_inline obj rational_plus( obj a, obj b )
{
  STUBBED_OUT(ZERO);
}

static _rs_inline obj rational_minus( obj a, obj b )
{
  STUBBED_OUT(ZERO);
}

static _rs_inline obj rational_mul( obj a, obj b )
{
  STUBBED_OUT(ZERO);
}

static _rs_inline obj rational_div( obj a, obj b )
{
  STUBBED_OUT(ZERO);
}


static int rational_cmp( obj a, obj b )
{
  STUBBED_OUT(0);
}

/************************ COMPLEX OPERATIONS ************************/

typedef struct _cmplx {
  obj re;
  obj im;
} cmplx;

static obj make_complex( obj re, obj im )
{
  return make2( rect_complex_class, re, im );
}

static _rs_inline cmplx real_to_complex( obj r )
{
  cmplx c;
  c.re = r;
  c.im = ZERO;
  return c;
}

static _rs_inline cmplx extract_complex( obj ch )
{
  cmplx c;
  c.re = gvec_ref(ch,SLOT(0));
  c.im = gvec_ref(ch,SLOT(1));
  return c;
}

#define RETURN_COMPLEX(re,im) if (EQ(im,ZERO)) return re; else return make_complex( re, im )

static _rs_inline obj complex_plus( cmplx a, cmplx b )
{
  obj re = basic_plus( a.re, b.re );
  obj im = basic_plus( a.im, b.im );

  RETURN_COMPLEX(re,im);
}

static _rs_inline obj complex_minus( cmplx a, cmplx b )
{
  obj re = basic_minus( a.re, b.re );
  obj im = basic_minus( a.im, b.im );

  RETURN_COMPLEX(re,im);
}

/*  (a + ib) * (c + id) = (ac - bd) + (ad + cb)i  */

static obj complex_mul( cmplx a_r, cmplx b_r )
{
  obj a = a_r.re;
  obj b = a_r.im;
  obj c = b_r.re;
  obj d = b_r.im;

  obj re = basic_minus( basic_mul( a, c ), basic_mul( b, d ) );
  obj im = basic_plus( basic_mul( a, d ), basic_mul( c, b ) );

  RETURN_COMPLEX(re,im);
}

/*  ??? */

static _rs_inline obj complex_div( cmplx a, cmplx b )
{
  obj re = basic_div( a.re, b.re );
  obj im = basic_div( a.im, b.im );

  RETURN_COMPLEX(re,im);
}

static _rs_inline obj complex_abs( cmplx a )
{
  return basic_plus( basic_mul( a.re, a.re ), basic_mul( a.im, a.im ) );
}

static int complex_cmp( cmplx a, cmplx b )
{
  return basic_cmp( complex_abs(a), complex_abs(b) );
}

/********************************************************************/

#include "numtower.ci"

#endif /* FULL_NUMERIC_TOWER */