[Easycalc-commits] PPCport/core/mlib bessels.cpp, NONE, 1.1 cmatrix.cpp, NONE, 1.1 cmatrix.h, NONE, 1.1 complex.cpp, NONE, 1.1 complex.h, NONE, 1.1 elliptic.cpp, NONE, 1.1 fl_num.cpp, NONE, 1.1 fl_num.h, NONE, 1.1 integ.cpp, NONE, 1.1 integ.h, NONE, 1.1 mathem.cpp, NONE, 1.1 mathem.h, NONE, 1.1 matrix.cpp, NONE, 1.1 matrix.h, NONE, 1.1 meqstack.cpp, NONE, 1.1 meqstack.h, NONE, 1.1 slist.cpp, NONE, 1.1 slist.h, NONE, 1.1 specfun.cpp, NONE, 1.1 specfun.h, NONE, 1.1 calcDB.cpp, 1.4, 1.5 calcDB.h, 1.1, 1.2 defuns.h, 1.1, 1.2 display.cpp, 1.4, 1.5 display.h, 1.3, 1.4 fp.cpp, 1.3, 1.4 fp.h, 1.1, 1.2 funcs.cpp, 1.4, 1.5 funcs.h, 1.3, 1.4 guess.cpp, 1.1, 1.2 history.cpp, 1.2, 1.3 history.h, 1.2, 1.3 konvert.cpp, 1.4, 1.5 konvert.h, 1.2, 1.3 stack.cpp, 1.4, 1.5 stack.h, 1.2, 1.3

[Mapi B. <ma...@us...>]

Update of /cvsroot/easycalc/PPCport/core/mlib
In directory 23jxhf1.ch3.sourceforge.com:/tmp/cvs-serv31293

Modified Files:
	calcDB.cpp calcDB.h defuns.h display.cpp display.h fp.cpp fp.h 
	funcs.cpp funcs.h guess.cpp history.cpp history.h konvert.cpp 
	konvert.h stack.cpp stack.h 
Added Files:
	bessels.cpp cmatrix.cpp cmatrix.h complex.cpp complex.h 
	elliptic.cpp fl_num.cpp fl_num.h integ.cpp integ.h mathem.cpp 
	mathem.h matrix.cpp matrix.h meqstack.cpp meqstack.h slist.cpp 
	slist.h specfun.cpp specfun.h 
Log Message:
1.25a Compiling, linking, running, limited function set, with bugs

--- NEW FILE: bessels.cpp ---
/*
 *   $Id: bessels.cpp,v 1.1 2009/10/17 13:48:34 mapibid Exp $
 * 
 *   Scientific Calculator for Palms.
 *   Copyright (C) 1999,2000 Ondrej Palkovsky
 *   This file Copyright (C) 2000 Rafael R. Sevilla
 *
 *   This program is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation; either version 2 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program; if not, write to the Free Software
 *   Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *   Boston, MA  02110-1301, USA.
 * 
 *  You can contact Ondrej at 'on...@pe...'.
 *  You can contact me at 'di...@pa...'.
 *
 *  Portions of this code are from the Cephes Math Library, Release 2.8
 *  Copyright (C) 1984, 1987, 1989, 2000 by Stephen L. Moshier
 */
#include "StdAfx.h"

#include "compat/PalmOS.h"
#include "compat/MathLib.h"
#include "konvert.h"
#include "funcs.h"
#include "mathem.h"
#include "stack.h"
#include "core/prefs.h"
#include "display.h"
#include "specfun.h"

static double p0[5] = {
   1.0000000000e0,
  -0.1098628627e-2,
   0.2734510407e-4,
  -0.2073370639e-5,
   0.2093887211e-6
}
;
static double q0[5] = {
  -0.1562499995e-1,
   0.1430488765e-3,
  -0.6911147651e-5, 
   0.7621095161e-6,
   0.9349451520e-7
};

static double p1[5] = {
   1.0000000000e0,
   0.1831050000e-2,
  -0.3516396496e-4,
   0.2457520174e-5,
  -0.2403370190e-6
};

static double q1[5] = {
   0.04687499995e0, 
  -0.20026908730e-3,
   0.84491990960e-5,
  -0.88228987000e-6,
   0.10578741200e-6
};

#define TWOOVRPI 0.63661977236758134307
#define PIOVR4   0.78539816339744830962
#define THRPIOV4 2.35619449019234492885

/***********************************************************************
 *
 * FUNCTION:     besselj01
 * 
 * DESCRIPTION:  Compute Bessel function of the first kind, order 0 or 1
 *
 * PARAMETERS:   n (either 0 or 1), x
 *
 * RETURN:       J_0(x) or J_1(x)
 *      
 ***********************************************************************/
static double besselj01(int n, double x) SPECFUN;
static double besselj01(int n, double x)
{
  double y, fact, pifact, ax, z, xx, rv;
  double *p, *q, *r, *s;
  static double r0[6] = { 57568490574.0, -13362590354.0, 651619640.7,
			 -11214424.18, 77392.33017, -184.9052456 };
  static double s0[6] = { 57568490411.0, 1029532985.0, 9494680.718,
			 59272.64853, 267.8532712, 1.0 };
  static double r1[6] = { 72362614232.0, -7895059235.0, 242396853.1,
			  -2972611.439, 15704.48260, -30.16036606 };
  static double s1[6] = { 144725228442.0, 2300535178.0, 18583304.74,
			  99447.43394, 376.9991397, 1.0 };

  if (n == 0) {
    p = p0;
    q = q0;
    r = r0;
    s = s0;
    fact = 1.0;
    pifact = PIOVR4; /* pi/4 */
  } else {
    p = p1;
    q = q1;
    r = r1;
    s = s1;
    fact = x;
    pifact = THRPIOV4; /* 3*pi/4 */
  }

  if (fabs(x) < 8) {
    y = x*x;
    return(fact*polevl(y, r, 6)/polevl(y, s, 6));
  }
  ax = fabs(x);
  z = 8/ax;
  y = z*z;
  xx = ax - pifact;
  rv = sqrt(TWOOVRPI/ax)*
    (cos(xx)*polevl(y, p, 5) - z*sin(xx)*polevl(y, q, 5));
  rv = (x < 0 && n == 1) ? -rv : rv;
  return(rv);
}

/***********************************************************************
 *
 * FUNCTION:     bessely01
 * 
 * DESCRIPTION:  Compute Bessel function of the second kind, order 0 or 1
 *
 * PARAMETERS:   n (either 0 or 1), x
 *
 * RETURN:       Y_0(x) or Y_1(x)
 *      
 ***********************************************************************/
static double bessely01(int n, double x) SPECFUN;
static double bessely01(int n, double x)
{
  double y, fact, pifact, onovx, z, xx, rv;
  double *p, *q, *r, *s;
  int sdeg;
  static double r0[6] = { -2957821389.0, 7062834065.0, -512359803.6,
			  10879881.29, -86327.92757, 228.4622733 };
  static double s0[6] = { 40076544269.0, 745249964.8, 7189466.438,
			  47447.26470, 226.1030244, 1.0 };
  static double r1[6] = { -0.4900604943e13, 0.1275274390e13, 
			  -0.5153438139e11, 0.7349264551e9,
			  -0.4237922726e7, 0.8511937935e4 };
  static double s1[7] = { 0.2499580570e14, 0.4244419664e12,
			  0.3733650367e10, 0.2245904002e8,
			  0.1020426050e6, 0.3549632885e3, 1.0 };
  if (n == 0) {
    p = p0;
    q = q0;
    r = r0;
    s = s0;
    fact = 1.0;
    sdeg = 6;
    pifact = PIOVR4; /* pi/4 */
    onovx = 0;
  } else {
    p = p1;
    q = q1;
    r = r1;
    s = s1;
    fact = x;
    sdeg= 7;
    pifact = THRPIOV4; /* 3*pi/4 */
    onovx = 1/x;
  }

  if (x < 8) {
    y = x*x;
    return(fact*polevl(y, r, 6)/polevl(y, s, sdeg) +
	   TWOOVRPI*(besselj01(n, x)*log(x) - onovx));
  }
  z = 8/x;
  y = z*z;
  xx = x - pifact;
  rv = sqrt(TWOOVRPI/x)*
    (sin(xx)*polevl(y, p, 5) + z*cos(xx)*polevl(y, q, 5));
  return(rv);
}

/***********************************************************************
 *
 * FUNCTION:     besseljn
 * 
 * DESCRIPTION:  Compute Bessel function of
 *               the first kind, order n
 *
 * PARAMETERS:   n, x, j01
 *
 * RETURN:       J_n(x)
 *      
 ***********************************************************************/
static double besseljn(int n, double x,
		       double (*j01)(int n, double x)) SPECFUN;
static double besseljn(int n, double x, double (*j01)(int n, double x))
{
  double pkm2, pkm1, pk, xk, r, ans;
  int k, sign;

  if( n < 0 ) {
    n = -n;
    if( (n & 1) == 0 )	/* -1**n */
      sign = 1;
    else
      sign = -1;
  } else
    sign = 1;

  if( x < 0.0 ) {
    if( n & 1 )
      sign = -sign;
    x = -x;
  }

  if( n == 0 )
    return( sign * j01(0, x) );
  if( n == 1 )
    return( sign * j01(1, x) );
  if( n == 2 )
    return( sign * (2.0 * j01(1, x) / x  -  j01(0, x)) );

  if( x < EPSILON )
    return( 0.0 );

  k = 53;

  pk = 2 * (n + k);
  ans = pk;
  xk = x * x;

  do {
    pk -= 2.0;
    ans = pk - (xk/ans);
  } while( --k > 0 );
  ans = x/ans;

  /* backward recurrence */

  pk = 1.0;
  pkm1 = 1.0/ans;
  k = n-1;
  r = 2 * k;

  do {
    pkm2 = (pkm1 * r  -  pk * x) / x;
    pk = pkm1;
    pkm1 = pkm2;
    r -= 2.0;
  } while( --k > 0 );

  if( fabs(pk) > fabs(pkm1) )
    ans = j01(1, x)/pk;
  else
    ans = j01(0, x)/pkm1;
  return( sign * ans );
}

/***********************************************************************
 *
 * FUNCTION:     besselyn
 * 
 * DESCRIPTION:  Compute Bessel function of
 *               the second kind, order n
 *
 * PARAMETERS:   n, x, y01
 *
 * RETURN:       Y_n(x)
 *      
 ***********************************************************************/
static double besselyn(int n, double x,
		       double (*y01)(int n, double x)) SPECFUN;
static double besselyn(int n, double x, double (*y01)(int n, double x))
{
  double an, anm1, anm2, r;
  int k, sign;

  if( n < 0 ) {
    n = -n;
    if( (n & 1) == 0 )	/* -1**n */
      sign = 1;
    else
      sign = -1;
  } else
    sign = 1;

  if( n == 0 )
    return( sign * y01(0, x) );
  if( n == 1 )
    return( sign * y01(1, x) );

  /* forward recurrence on n */

  anm2 = y01(0, x);
  anm1 = y01(1, x);
  k = 1;
  r = 2 * k;
  do {
    an = r * anm1 / x  -  anm2;
    anm2 = anm1;
    anm1 = an;
	r += 2.0;
	++k;
  } while( k < n );

  return( sign * an );
}

/***********************************************************************
 *
 * FUNCTION:     besseli01
 * 
 * DESCRIPTION:  Compute modified Bessel function of the first kind,
 *               order 0 or 1
 *
 * PARAMETERS:   n (0 or 1 only), x
 *
 * RETURN:       I_0(x) or I_1(x)
 *      
 ***********************************************************************/
static double besseli01(int n, double x) SPECFUN;
static double besseli01(int n, double x)
{
  double y, ax, fact;
  double *p, *q;
  static double p0[7] = { 1.0, 3.5156229, 3.0899424, 1.2067492, 0.2659732,
			0.360768e-1, 0.45813e-2 };
  static double q0[9] = { 0.39894228, 0.1328592e-1, 0.225319e-2, -0.157565e-2,
			0.916281e-2, -0.2057706e-1, 0.2635537e-1,
			-0.1647633e-1, 0.392377e-2 };
  static double p1[7] = { 0.5, 0.87890594, 0.51498869, 0.15084934,
			  0.2658733e-1, 0.301532e-2, 0.32411e-3 };
  static double q1[9] = { 0.39894228, -0.3988024e-1, -0.362018e-2, 0.163801e-2,
			 -0.1031555e-1, 0.2282967e-1, -0.2895312e-1,
			 0.1787654e-1, -0.420059e-2 };

  p = (n == 0) ? p0 : p1;
  q = (n == 0) ? q0 : q1;

  if (fabs(x) < 3.75) {
    fact = (n == 0) ? 1 : x;
    y = (x*x/14.0625);
    return(fact*polevl(y, p, 7));
  }
  ax = fabs(x);
  y = 3.75/ax;
  return((exp(ax)/sqrt(ax))*polevl(y, q, 9));
}

/***********************************************************************
 *
 * FUNCTION:     besselk01
 * 
 * DESCRIPTION:  Compute modified Bessel function of the 2nd kind,
 *               order 0 or 1.
 *
 * PARAMETERS:   n (either 0 or 1 only), x
 *
 * RETURN:       K_0(x) or K_1(x)
 *      
 ***********************************************************************/
static double besselk01(int n, double x) SPECFUN;
static double besselk01(int n, double x)
{
  double y, fact;
  double *p, *q;
  double p0[7] = { -0.57721566, 0.42278420, 0.23069756, 0.3488590e-1,
		  0.262698e-2, 0.10750e-3, 0.74e-5 };
  double q0[7] = { 1.25331414, -0.7832358e-1, 0.2189568e-1, -0.1062446e-1,
		  0.587872e-2, -0.251540e-2, 0.53208e-3 };
  double p1[7] = { 1.0, 0.15443144, -0.67278579, -0.18156897, -0.1919402e-1,
		  -0.110404e-2, -0.4686e-4 };
  double q1[7] = { 1.25331414, 0.23498619, -0.3655620e-1, 0.1504268e-1,
		  -0.780353e-2, 0.325614e-2, -0.68245e-3 };
  int sign;

  if (n == 0) {
    p = p0;
    q = q0;
    sign = -1;
    fact = 1.0;
  } else {
    p = p1;
    q = q1;
    sign = 1;
    fact = 1.0/x;
  }
  if (x <= 2.0) {
    y = x*x/4.0;
    return((sign*log(x/2.0)*besseli01(n, x))+
	   fact*polevl(y, p, 7));
  }
  y = 2.0/x;
  return((exp(-x)/sqrt(x))*polevl(y, q, 7));
}

/***********************************************************************
 *
 * FUNCTION:     besselin
 * 
 * DESCRIPTION:  Compute modified Bessel function of the first kind,
 *               order n
 *
 * PARAMETERS:   n, x
 *
 * RETURN:       I_n(x)
 *      
 ***********************************************************************/
#define IACC 40
#define BIGNO 1e10
#define BIGNI 1e-10
static double besselin(int n, double x,
		       double (*i01)(int n, double x)) SPECFUN;
static double besselin(int n, double x, double (*i01)(int n, double x))
{
  double tox, bip, bi, bim, retval;
  int j, m;

  if (n < 0)
    n = -n;
  if (n == 0 || n == 1)
    return(i01(n, x));
  if (x == 0)
    return(0);
  tox = 2.0/fabs(x);
  bip = 0.0;
  bi = 1.0;
  retval = 0.0;
  m = 2*(n+(int)(sqrt(IACC*n)));
  for (j=m; j>=1; j--) {
    bim = bip + j*tox*bi;
    bip = bi;
    bi = bim;
    if (fabs(bi) > BIGNO) {
      retval *= BIGNI;
      bi *= BIGNI;
      bip *= BIGNI;
    }
    if (j == n)
      retval = bip;
  }
  retval *= i01(0, x)/bi;
  if (x < 0 && (n & 1))
    return(-retval);
  return(retval);
}

/***********************************************************************
 *
 * FUNCTION:     besselkn
 * 
 * DESCRIPTION:  Compute modified Bessel function of the second kind,
 *               order n
 *
 * PARAMETERS:   n, x
 *
 * RETURN:       K_n(x)
 *      
 ***********************************************************************/
static double besselkn(int n, double x,
		       double (*k01)(int n, double x)) SPECFUN;
static double besselkn(int n, double x, double (*k01)(int n, double x))
{
  double tox, bkm, bkp, bk;
  int j;

  if (n < 0)
    n = -n;

  if (n == 0 || n == 1)
    return(k01(n, x));
  tox = 2.0/x;
  bkm = k01(0, x);
  bk = k01(1, x);
  for (j=1; j<n; j++) {
    bkp = bkm + j*tox*bk;
    bkm = bk;
    bk = bkp;
  }
  return(bk);
}

/***********************************************************************
 *
 * FUNCTION:     math_bessels
 * 
 * DESCRIPTION:  Compute Bessel and related functions
 *
 * PARAMETERS:   On stack - 2 numbers
 *
 * RETURN:       On stack - 1 number
 *      
 ***********************************************************************/
CError
math_bessels(Functype *func, CodeStack *stack)
{

  double x, n, res;
  CError err;

  err = stack_get_val(stack, &x, real);
  if (err)
    return(err);

  switch (func->num) {
  case MATH_BESSELJ:
  case MATH_BESSELY:
    err = stack_get_val(stack, &n, real);
    if (err)
      return(err);
    /* integer bessels only */
    if (!IS_ZERO(n - round(n)) ||
	(IS_ZERO(x) && (func->num == MATH_BESSELY)))
      return(c_badarg);
    res = (func->num == MATH_BESSELJ) ? besseljn((int) n, x, besselj01) : 
      besselyn((int) n, x, bessely01);
    break;
  case MATH_BESSELI:
  case MATH_BESSELK:
    err = stack_get_val(stack, &n, real);
    if (err)
      return(err);
    /* integer bessels only */
    if (!IS_ZERO(n - round(n)) ||
	(IS_ZERO(x) && (func->num == MATH_BESSELK)))
      return(c_badarg);
    res = (func->num == MATH_BESSELI) ? besselin((int) n, x, besseli01) :
      besselkn((int) n, x, besselk01);
    break;
  case MATH_AIRY:
  case MATH_BIRY:
  default:
    return(c_internal);
  }
  return stack_add_val(stack, &res, real);
}


--- NEW FILE: fl_num.h ---
/*
 *   Scientific Calculator for Palms.
 *   Copyright (C) 2000 Ondrej Palkovsky
 *
 *   This program is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation; either version 2 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program; if not, write to the Free Software
 *   Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *   Boston, MA  02110-1301, USA.
 * 
 *  You can contact me at 'on...@pe...'.
*/

Int16 fl_num(TCHAR *input,double *cislo) MLIB;

Index: guess.cpp
===================================================================
RCS file: /cvsroot/easycalc/PPCport/core/mlib/guess.cpp,v
retrieving revision 1.1
retrieving revision 1.2
diff -C2 -d -r1.1 -r1.2
*** guess.cpp	24 Sep 2009 21:16:50 -0000	1.1
--- guess.cpp	17 Oct 2009 13:48:34 -0000	1.2
***************
*** 28,32 ****
  //#include <StringMgr.h>
  
! #include "segment.h"
  #include "compat/MathLib.h"
  #include "stack.h"
--- 28,32 ----
  //#include <StringMgr.h>
  
! #include "compat/segment.h"
  #include "compat/MathLib.h"
  #include "stack.h"
***************
*** 454,458 ****
  
      if (item.type == real) {
!         result = (TCHAR *) MemPtrNew(MAX_FP_NUMBER+1);
          guessed = guess_real(item.u.realval,result);
          if (!guessed) {
--- 454,458 ----
  
      if (item.type == real) {
!         result = (TCHAR *) MemPtrNew((MAX_FP_NUMBER+1)*sizeof(TCHAR));
          guessed = guess_real(item.u.realval,result);
          if (!guessed) {
***************
*** 463,467 ****
      }
      else if (item.type == complex) {
!         result = (TCHAR *) MemPtrNew(MAX_FP_NUMBER*3);
          if (IS_ZERO(item.u.cplxval->imag)){
             if (guess_real(item.u.cplxval->real,result))
--- 463,467 ----
      }
      else if (item.type == complex) {
!         result = (TCHAR *) MemPtrNew((MAX_FP_NUMBER*3)*sizeof(TCHAR));
          if (IS_ZERO(item.u.cplxval->imag)){
             if (guess_real(item.u.cplxval->real,result))


--- NEW FILE: specfun.h ---
/*   
 *   $Id: specfun.h,v 1.1 2009/10/17 13:48:34 mapibid Exp $
 * 
 *   Scientific Calculator for Palms.
 *   This file Copyright (C) 2000 Rafael Sevilla
 *
 *   This program is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation; either version 2 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program; if not, write to the Free Software
 *   Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *   Boston, MA  02110-1301, USA.
 * 
 *  You can contact me at 'di...@pa...'.
 *
 */

#ifndef _SPECFUN_H
#define _SPECFUN_H

#include "compat/segment.h"

CError math_gambeta(Functype *func, CodeStack *stack) SPECFUN;
CError math_bessels(Functype *func, CodeStack *stack) SPECFUN;
CError math_elliptic(Functype *func, CodeStack *stack) SPECFUN;
CError math_hypergeo(Functype *func, CodeStack *stack) SPECFUN;
CError math_orthpol(Functype *func, CodeStack *stack) SPECFUN;
double polevl(double x, double *coef, int N ) SPECFUN;

#define EPSILON 1e-10

#endif


Index: calcDB.cpp
===================================================================
RCS file: /cvsroot/easycalc/PPCport/core/mlib/calcDB.cpp,v
retrieving revision 1.4
retrieving revision 1.5
diff -C2 -d -r1.4 -r1.5
*** calcDB.cpp	24 Sep 2009 21:16:50 -0000	1.4
--- calcDB.cpp	17 Oct 2009 13:48:34 -0000	1.5
***************
*** 22,27 ****
   *  You can contact me at 'on...@pe...'.
  */
- 
- 
  #include "stdafx.h"
  
--- 22,25 ----
***************
*** 100,104 ****
                  err = c_memory;
              DmReleaseRecord(gDB,recordNumber-1,true);
!             return err;
          }
      }
--- 98,102 ----
                  err = c_memory;
              DmReleaseRecord(gDB,recordNumber-1,true);
!             return (err);
          }
      }
***************
*** 109,117 ****
          MemHandleUnlock(myRecordMemHandle);
          DmReleaseRecord(gDB,recordNumber,true); /* Now the recordNumber contains real index */
!     }
!     else
          err = c_memory;
  
!     return err;
  }
  
--- 107,114 ----
          MemHandleUnlock(myRecordMemHandle);
          DmReleaseRecord(gDB,recordNumber,true); /* Now the recordNumber contains real index */
!     } else
          err = c_memory;
  
!     return (err);
  }
  
***************
*** 184,192 ****
      dbPackedStack *pstack;
      Int16 psize,rsize;
!     TCHAR *ptr;
      CError err;
  
      pstack = stack_pack(stack,&psize);
!     record = (dbPackedRecord *) MemPtrNew(sizeof(*record)+psize+StrLen(origtext)+1);
      if (!record)
          return c_memory;
--- 181,190 ----
      dbPackedStack *pstack;
      Int16 psize,rsize;
!     byte *ptr;
      CError err;
  
      pstack = stack_pack(stack,&psize);
!     rsize = sizeof(*record) + psize + (StrLen(origtext)+1) * sizeof(TCHAR);
!     record = (dbPackedRecord *) MemPtrNew(rsize);
      if (!record)
          return c_memory;
***************
*** 198,206 ****
      StrCopy(record->u.stack.header.param_name,parameter_name);
  
! //    ptr = &record->u.stack.data;
!     ptr = record->u.stack.data;
      ptr += record->u.stack.header.orig_offset;
!     StrCopy(ptr,origtext);
!     rsize = sizeof(*record)+psize+StrLen(origtext)+1;
      rsize -= sizeof(dbStackItem)>sizeof(dbPackedStack)?sizeof(dbStackItem):sizeof(dbPackedStack);
      err = db_create_record(record,rsize);
--- 196,203 ----
      StrCopy(record->u.stack.header.param_name,parameter_name);
  
! //    ptr = (byte *) &(record->u.stack.data);
!     ptr = (byte *) (record->u.stack.data);
      ptr += record->u.stack.header.orig_offset;
!     StrCopy((TCHAR *) ptr, origtext);
      rsize -= sizeof(dbStackItem)>sizeof(dbPackedStack)?sizeof(dbStackItem):sizeof(dbPackedStack);
      err = db_create_record(record,rsize);
***************
*** 377,381 ****
               && record->u.data[0].rpn.type == type)) {
              result->type[count] = record->type;
!             result->list[count] = (TCHAR *) MemPtrNew(StrLen(record->name)+1);
              StrCopy(result->list[count],record->name);
              count++;
--- 374,378 ----
               && record->u.data[0].rpn.type == type)) {
              result->type[count] = record->type;
!             result->list[count] = (TCHAR *) MemPtrNew((StrLen(record->name)+1)*sizeof(TCHAR));
              StrCopy(result->list[count],record->name);
              count++;
***************
*** 413,420 ****
          return c_badarg;
      }
!     *dest = (TCHAR *) MemPtrNew(StrLen(&tmp->u.stack.data[tmp->u.stack.header.orig_offset])+1);
!     StrCopy(*dest,&tmp->u.stack.data[tmp->u.stack.header.orig_offset]);
      if (param)
!       StrCopy(param,tmp->u.stack.header.param_name);
  
      MemPtrFree(tmp);
--- 410,420 ----
          return c_badarg;
      }
!     *dest = (TCHAR *) MemPtrNew((StrLen((TCHAR *) (tmp->u.stack.data + tmp->u.stack.header.orig_offset))
!                                  +1
!                                 )*sizeof(TCHAR)
!                                );
!     StrCopy(*dest, (TCHAR *) (tmp->u.stack.data + tmp->u.stack.header.orig_offset));
      if (param)
!         StrCopy(param, tmp->u.stack.header.param_name);
  
      MemPtrFree(tmp);

--- NEW FILE: integ.cpp ---
/*
 *   $Id: integ.cpp,v 1.1 2009/10/17 13:48:34 mapibid Exp $
 * 
 *   Scientific Calculator for Palms.
 *   Copyright (C) 1999,2000 Ondrej Palkovsky
 *
 *   This program is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation; either version 2 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program; if not, write to the Free Software
 *   Foundation, Inc., 51 Franklin Street, Fifth Floor,
[...989 lines suppressed...]
		stack_delete(f1);
		stack_delete(f2);
		goto error;
	}

	result = integ_intersect(min,max,f1,f2,error,argarr);
	stack_delete(f1);
	stack_delete(f2);
	
	if (!finite(result)) {
		err = c_compimp;
		goto error;
	}
	err = stack_add_val(stack,&result,real);
error:
	if (argarr)
		stack_delete(argarr);
	return err;
}


Index: konvert.h
===================================================================
RCS file: /cvsroot/easycalc/PPCport/core/mlib/konvert.h,v
retrieving revision 1.2
retrieving revision 1.3
diff -C2 -d -r1.2 -r1.3
*** konvert.h	24 Jun 2009 21:39:59 -0000	1.2
--- konvert.h	17 Oct 2009 13:48:34 -0000	1.3
***************
*** 27,31 ****
  #define _konvert_h_
  
! #include "segment.h"
  
  #define MAX_FUNCNAME 10
--- 27,31 ----
  #define _konvert_h_
  
! #include "compat/segment.h"
  
  #define MAX_FUNCNAME 10
***************
*** 33,41 ****
  
  /* First function/variable character */
! #define IS_FUNC_1LETTER(x) (((x)>='a' && (x)<='z') || (x)=='_' || (x)==(TCHAR)181 || \
!                            ((x)>='G' && (x)<='Z'))
  /* Other characters, inside varname can be numbers */
! #define IS_FUNC_LETTER(x)  (IS_FUNC_1LETTER(x) || ((x)>='0' && (x)<='9') || \
! 								((x)>='A' && (x)<='Z'))
  
  typedef enum {
--- 33,41 ----
  
  /* First function/variable character */
! #define IS_FUNC_1LETTER(x) (((x)>=_T('a') && (x)<=_T('z')) || (x)==_T('_') || (x)==(TCHAR)181 || \
!                            ((x)>=_T('G') && (x)<=_T('Z')))
  /* Other characters, inside varname can be numbers */
! #define IS_FUNC_LETTER(x)  (IS_FUNC_1LETTER(x) || ((x)>=_T('0') && (x)<=_T('9')) || \
! 								((x)>=_T('A') && (x)<=_T('Z')))
  
  typedef enum {
***************
*** 116,120 ****
  typedef struct {
  	UInt16 size;
! 	Complex *item;
  }List;
  
--- 116,120 ----
  typedef struct {
  	UInt16 size;
! 	Complex item[1];
  }List;
  
***************
*** 122,131 ****
  typedef struct {
  	Int16 rows,cols;
! 	double *item;
  }Matrix;
  
  typedef struct {
  	Int16 rows,cols;
! 	Complex *item;
  }CMatrix;
  	
--- 122,131 ----
  typedef struct {
  	Int16 rows,cols;
! 	double item[1];
  }Matrix;
  
  typedef struct {
  	Int16 rows,cols;
! 	Complex item[1];
  }CMatrix;
  	

Index: stack.h
===================================================================
RCS file: /cvsroot/easycalc/PPCport/core/mlib/stack.h,v
retrieving revision 1.2
retrieving revision 1.3
diff -C2 -d -r1.2 -r1.3
*** stack.h	24 Jun 2009 21:39:59 -0000	1.2
--- stack.h	17 Oct 2009 13:48:34 -0000	1.3
***************
*** 27,31 ****
  
  #include "calcDB.h"
! #include "segment.h"
  
  dbPackedStack * stack_pack(CodeStack *stack,Int16 *size) MLIB;
--- 27,31 ----
  
  #include "calcDB.h"
! #include "compat/segment.h"
  
  dbPackedStack * stack_pack(CodeStack *stack,Int16 *size) MLIB;

Index: calcDB.h
===================================================================
RCS file: /cvsroot/easycalc/PPCport/core/mlib/calcDB.h,v
retrieving revision 1.1
retrieving revision 1.2
diff -C2 -d -r1.1 -r1.2
*** calcDB.h	24 Jun 2009 21:44:22 -0000	1.1
--- calcDB.h	17 Oct 2009 13:48:34 -0000	1.2
***************
*** 32,36 ****
      Int16 size;
      rpntype *type;
!     TCHAR **list;
      /* Do not add anything here, it gets overwritten */
  }dbList;
--- 32,36 ----
      Int16 size;
      rpntype *type;
!     TCHAR *list[1];
      /* Do not add anything here, it gets overwritten */
  }dbList;
***************
*** 45,54 ****
      Trpn rpn;
      UInt16 datasize; /* size of additional packed data */
!     TCHAR *data;
  }dbStackItem;
  
  typedef struct {
      dbStackHeader header;
!     TCHAR *data;
  }dbPackedStack;
  
--- 45,54 ----
      Trpn rpn;
      UInt16 datasize; /* size of additional packed data */
!     byte data[1];
  }dbStackItem;
  
  typedef struct {
      dbStackHeader header;
!     byte data[1];
  }dbPackedStack;
  

--- NEW FILE: fl_num.cpp ---
/*   
 *   $Id: fl_num.cpp,v 1.1 2009/10/17 13:48:34 mapibid Exp $
 * 
 *   Scientific Calculator for Palms.
 *   Copyright (C) 1999 Ondrej Palkovsky
 *
 *   This program is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation; either version 2 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program; if not, write to the Free Software
 *   Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *   Boston, MA  02110-1301, USA.
 * 
 *  You can contact me at 'on...@pe...'.
*/
#include "StdAfx.h"

#include "compat/PalmOS.h"

#include "compat/segment.h"
#include "konvert.h"
#include "display.h"
#include "fp.h"
#include "core/prefs.h"

#ifdef PALMOS
#define isdigit(x) ((x>=_T('0')) && (x<=_T('9')))
#endif

/* Is a number, equal to /[0-9]+/ */
static Int16 is_ccbz(TCHAR *input,double *cislo) MLIB;
static Int16
is_ccbz(TCHAR *input,double *cislo)
{
	Int16 i,j;
	Int16 spaces = 0;

	/* Count to 'i' the length of a number */
	for (i=0;;i++) {
		if (input[i]>=_T('0') && input[i]<=_T('9'))
			j = input[i]-_T('0');
		else if (input[i]>=_T('A') && input[i]<=_T('F'))
			j = input[i] - _T('A') + 10;
		else if (calcPrefs.acceptOSPref && input[i] == flSpaceChar) {
			spaces++;
			continue;
		} else
			break;
		if (j>dispPrefs.base-1)
			break;
	}		

	if (cislo) {
		*cislo=0.0;
		for (j=0;j<i;j++) {
			if (input[j] == flSpaceChar)
				continue;
			*cislo*=dispPrefs.base;
			*cislo+=input[j]<=_T('9')?input[j]-_T('0'):input[j]-_T('A')+10;
		}
		j= (Int16) (*cislo);
	}

	if (spaces == i)
		return 0;
	return i;
}

static Int16 has_expon(TCHAR *input,double *cislo) MLIB;
static Int16
has_expon(TCHAR *input,double *cislo)
{
	TCHAR *tmp=input;
	Int16 delka=0;
	Int16 i;
	
	if (*tmp==_T('E')) {
		Int16 negate=0;
		double expon;
		
		tmp++;delka++;
		if (*tmp==_T('-')) {
			negate=1;
			tmp++;delka++;
		}
		else if (*tmp==_T('+')) {
			tmp++;delka++;
		}
	        if (!(i=is_ccbz(tmp,&expon)))
		  return 0;
		delka+=i;
/* Prevent excessive waiting */
		if (expon>350.0)
		  return 0;
		
		*cislo=1.0;
		if (!negate)
		  for (i=0;i<expon;i++)
		    *cislo*=10.0;
		else
		  for (i=0;i<expon;i++)
		    *cislo/=10.0;
	}
	return delka;
}

/* If floating number, return number of characters describing
 * the number */
Int16 fl_num(TCHAR *input,double *cislo) MLIB;
Int16
fl_num(TCHAR *input,double *cislo)
{	
	TCHAR *tmp=input;
	double cele;
	double expon;
	Int16 i=0,j;
	double dfl;
	Int16 delka=0;
		
	if ((i=is_ccbz(tmp,&cele))) {
		tmp+=i;
		delka+=i;
	}
	else
	  cele=0;       
	
	if (calcPrefs.acceptOSPref && *tmp != flPointChar && delka==0)
		return 0;
	else if (!calcPrefs.acceptOSPref && *tmp != _T('.') && *tmp!=_T(',') && delka==0)
		return 0;
	else if (((calcPrefs.acceptOSPref && *tmp != flPointChar)
		  || (!calcPrefs.acceptOSPref && *tmp!=_T('.') && *tmp!=_T(',')))
		 && delka>0) {
		*cislo=cele;
		if ((i=has_expon(tmp,&expon))) {
			*cislo*=expon;
			delka+=i;tmp+=i;	
		}

		/* If a space is a thousend delimiter, check that
		 * there is not another number following a few spaces.
		 * this would interpret '3 4' as '34' or '3*4' depending
		 * on preferences if I didn't check here
		 */
		if (flSpaceChar == _T(' ') && *tmp == _T(' ')) {
			for (;*tmp == ' ';tmp++)
				;
			if (*tmp >= _T('0') && *tmp <= _T('9'))
				/* space and number found */
				return 0;
		}
		return delka;
	}
	tmp++;delka++;
	
	/* If we have floating point part */
	if ((i=is_ccbz(tmp,NULL))) {
		delka+=i;
		
		dfl=1.0/(double)dispPrefs.base;
		for (j=0;j<i;j++) {		
			if (calcPrefs.acceptOSPref && tmp[j] == flSpaceChar)
				continue;
			cele+=((double)(tmp[j]<=_T('9')?tmp[j]-_T('0'):tmp[j]-_T('A')+10))*dfl;
			dfl/=dispPrefs.base;
		}
		tmp+=i;
	}
	if (!delka)
		return 0;
	*cislo=cele;

	if ((i=has_expon(tmp,&expon))) {
		*cislo*=expon;
		delka+=i;tmp+=i;	
	}

	/* If a space is a thousend delimiter, check that
	 * there is not another number following a few spaces.
	 * this would interpret '3 4' as '34' or '3*4' depending
	 * on preferences if I didn't check here
	 */
	if (flSpaceChar == _T(' ') && *tmp == _T(' ')) {
		for (;*tmp == _T(' ');tmp++)
			;
		if (*tmp >= _T('0') && *tmp <= _T('9'))
			/* space and number found */
			return 0;
	}
	return delka;
}

--- NEW FILE: integ.h ---
/*   
 *   $Id: integ.h,v 1.1 2009/10/17 13:48:34 mapibid Exp $
 * 
 *   Scientific Calculator for Palms.
 *   Copyright (C) 1999,2000 Ondrej Palkovsky
 *
 *   This program is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation; either version 2 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program; if not, write to the Free Software
 *   Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *   Boston, MA  02110-1301, USA.
 * 
 *  You can contact me at 'on...@pe...'.
*/

#ifndef _INTEG_H_
#define _INTEG_H_

/* Precision for solving functions */
#define DEFAULT_ERROR (1E-8)
#define DEFAULT_ROMBERG  6
#define MAX_ITER 2000

CError integ_fsolve(Functype *func,CodeStack *stack) NEWFUNC;
CError integ_fsimps(Functype *func,CodeStack *stack) NEWFUNC;
CError integ_fdydx(Functype *func,CodeStack *stack) NEWFUNC;
double integ_simps(double min, double max,CodeStack *fcstack,
		   double error,CodeStack *argarr) NEWFUNC;
double integ_zero(double min, double max, double value, CodeStack *fcstack, 
		  double error,Int16 funcnum,CodeStack *argarr) NEWFUNC;
double integ_derive1(double x,CodeStack *fcstack,double error,
		     CodeStack *argarr) NEWFUNC;
double integ_derive2(double x,CodeStack *fcstack,double error,
		     CodeStack *argarr) NEWFUNC;
double integ_romberg(double min,double max,CodeStack *fcstack,Int16 n,
		     CodeStack *argarr) NEWFUNC;
CError integ_fromberg(Functype *func,CodeStack *stack) NEWFUNC;
double integ_intersect(double min,double max, CodeStack *f1, CodeStack *f2,
		       double error,CodeStack *argarr) NEWFUNC;
CError integ_fintersect(Functype *func,CodeStack *stack) NEWFUNC;

#endif

Index: funcs.cpp
===================================================================
RCS file: /cvsroot/easycalc/PPCport/core/mlib/funcs.cpp,v
retrieving revision 1.4
retrieving revision 1.5
diff -C2 -d -r1.4 -r1.5
*** funcs.cpp	24 Sep 2009 21:16:50 -0000	1.4
--- funcs.cpp	17 Oct 2009 13:48:34 -0000	1.5
***************
*** 53,82 ****
  
  struct {
!         Trpn *argument;
!         Int16 argcount;
  }func_arguments[MAX_RECURS*2];
  Int16 func_argcount;
  
  const TCHAR *funcInterDef[] = {
!         /* covariance */
!         _T("sum((x(1)-mean(x(1)))*(x(2)-mean(x(2))))/dim(x(1))"),
[...1413 lines suppressed...]
!         item = stack_pop(stack);
!         ritem.allocsize = 0;
!         ritem.type = litem;
!         StrCopy(ritem.u.litemval.name,item.u.varname);
!         ritem.u.litemval.row = arg1;
!         ritem.u.litemval.col = arg2;
!         stack_push(stack,ritem);
!         rpn_delete(item);
!     } else switch (type2){
!         case list:
!             return list_item(func,stack);
!         case matrix:
!             return matrix_item(func,stack);
!         case cmatrix:
!             return cmatrix_item(func,stack);
!         default:
!             return c_badarg;
!     }
      return c_noerror;
  }

--- NEW FILE: cmatrix.h ---
/*
 *   $Id: cmatrix.h,v 1.1 2009/10/17 13:48:34 mapibid Exp $
 *
 *   Scientific Calculator for Palms.
 *   Copyright (C) 1999,2000,2001,2002 Ondrej Palkovsky
 *
 *   This program is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation; either version 2 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program; if not, write to the Free Software
 *   Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *   Boston, MA  02110-1301, USA.
 * 
 *  You can contact me at 'on...@pe...'.
*/

#ifndef _CMATRIX_H_
#define _CMATRIX_H_

#define CMATRIX_INPUT   2101
#define CMATRIX_DET     2102
#define CMATRIX_ECHELON 2103
#define CMATRIX_IDENTITY 2104
#define CMATRIX_QRQ 2105
#define CMATRIX_QRR 2106
#define CMATRIX_QRS 2107

CError cmatrix_input(Functype *func,CodeStack *stack) MATFUNC;
CMatrix * cmatrix_new(Int16 rows,Int16 cols) MATFUNC;
CMatrix * cmatrix_dup(CMatrix *m) MATFUNC;
void cmatrix_delete(CMatrix *m) MATFUNC;
Matrix * cmatrix_to_matrix(CMatrix *cm) MATFUNC;
CError cmatrix_oper(Functype *func,CodeStack *stack) MATFUNC;
CMatrix * matrix_to_cmatrix(Matrix *m) MATFUNC;
CError cmatrix_oper_cplx(Functype *func,CodeStack *stack) MATFUNC;
CError cmatrix_func2(Functype *func,CodeStack *stack) MATFUNC;
CError cmatrix_func(Functype *func,CodeStack *stack) MATFUNC;
CError cmatrix_dim(Functype *func, CodeStack *stack) MATFUNC;
CError cmatrix_item(Functype *func, CodeStack *stack) MATFUNC;

#endif

--- NEW FILE: mathem.cpp ---
/*
 *   $Id: mathem.cpp,v 1.1 2009/10/17 13:48:34 mapibid Exp $
 * 
 *   Scientific Calculator for Palms.
 *   Copyright (C) 1999,2000 Ondrej Palkovsky
 *
 *   This program is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation; either version 2 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program; if not, write to the Free Software
 *   Foundation, Inc., 51 Franklin Street, Fifth Floor,
[...1370 lines suppressed...]
		for(i=2; i < 65536 && i*i <= val; i++) {
			if (val % i == 0) {
				isprime = false;
				break;
			}
		}
		if (func->num == MATH_ISPRIME) {
			val = isprime ? 1 : 0;
			break;
		} else if (isprime) {
			break;
		} else if (func->num == MATH_NEXTPRIME) {
			val++;
		} else {
			val--;
		}
	}

	return stack_add_val(stack, &val, integer);
}

Index: stack.cpp
===================================================================
RCS file: /cvsroot/easycalc/PPCport/core/mlib/stack.cpp,v
retrieving revision 1.4
retrieving revision 1.5
diff -C2 -d -r1.4 -r1.5
*** stack.cpp	24 Sep 2009 21:16:50 -0000	1.4
--- stack.cpp	17 Oct 2009 13:48:34 -0000	1.5
***************
*** 321,325 ****
  			return c_badarg;
  		}
! 		tmp = (TCHAR *) MemPtrNew(StrLen(tmpitem.u.stringval)+1);
  		StrCopy(tmp,tmpitem.u.stringval);
  		*((TCHAR **)arg1)=tmp;
--- 321,325 ----
  			return c_badarg;
  		}
! 		tmp = (TCHAR *) MemPtrNew((StrLen(tmpitem.u.stringval)+1)*sizeof(TCHAR));
  		StrCopy(tmp,tmpitem.u.stringval);
  		*((TCHAR **)arg1)=tmp;
***************
*** 479,483 ****
  		}
  	} else if (reqtype == string) {
! 		value.allocsize = StrLen((TCHAR *)arg)+1;
  		value.u.stringval = (TCHAR *) MemPtrNew(value.allocsize);
  		StrCopy(value.u.stringval,(TCHAR *)arg);
--- 479,483 ----
  		}
  	} else if (reqtype == string) {
! 		value.allocsize = (StrLen((TCHAR *)arg)+1)*sizeof(TCHAR);
  		value.u.stringval = (TCHAR *) MemPtrNew(value.allocsize);
  		StrCopy(value.u.stringval,(TCHAR *)arg);

--- NEW FILE: meqstack.h ---
/*
 *   $Id: meqstack.h,v 1.1 2009/10/17 13:48:34 mapibid Exp $
 * 
 *   Scientific Calculator for Palms.
 *   Copyright (C) 2000 Ondrej Palkovsky
 *
 *   This program is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation; either version 2 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program; if not, write to the Free Software
 *   Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *   Boston, MA  02110-1301, USA.
 * 
 *  You can contact me at 'on...@pe...'.
*/

#ifndef _MEQSTACK_H_
#define _MEQSTACK_H_

#include "konvert.h"

typedef struct {
	Int16 allocated;
	Int16 count;
	Tmeq *array;
}Tdynmeq;

Int16 meq_push(Tdynmeq *array,Tmeq item);
void meq_free(Tdynmeq *array);
void meq_update_t_pcount(Tdynmeq *array,Int16 position,Int16 pcount);
void meq_update_f_pcount(Tdynmeq *array,Int16 position,Int16 pcount);
Tmeq meq_pop(Tdynmeq *array);
Tmeq meq_last(Tdynmeq *array);
Tmeq meq_fetch(Tdynmeq *array);
Int16 meq_count(Tdynmeq *array);
Tdynmeq *meq_new(void);


#endif

--- NEW FILE: specfun.cpp ---
/*
 *   $Id: specfun.cpp,v 1.1 2009/10/17 13:48:34 mapibid Exp $
 * 
 *   Scientific Calculator for Palms.
 *   Copyright (C) 1999,2000 Ondrej Palkovsky
 *   This file Copyright (C) 2000 Rafael R. Sevilla
 *
 *   This program is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation; either version 2 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program; if not, write to the Free Software
 *   Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *   Boston, MA  02110-1301, USA.
 * 
 *  You can contact Ondrej at 'on...@pe...'.
 *  You can contact me at 'di...@pa...'.
 *
 */
#include "StdAfx.h"

#include "compat/PalmOS.h"
#include "compat/MathLib.h"
#include "konvert.h"
#include "funcs.h"
#include "mathem.h"
#include "stack.h"
#include "core/prefs.h"
#include "display.h"
#include "specfun.h"

extern double euler_gamma(double z);

/* Arrays are indexed from 1, that's why the 0.0 in the beginning */
/* Coefficients of rational functions for ln(gamma) approxiamtion*/
static const double r1[] = {0.0, -2.66685511495E0, -2.44387534237E1,
		      -2.19698958928E1, 1.11667541262E1,
		      3.13060547623E0, 6.07771387771E-1,
		      1.19400905721E1, 3.14690115749E1,
		      1.52346874070E1};

static const double r2[] = {0.0, -7.83359299449E1, -1.42046296688E2,
		      1.37519416416E2, 7.86994924154E1,
		      4.16438922228E0, 4.70668766060E1,
		      3.13399215894E2,  2.63505074721E2,
		      4.33400022514E1};

static const double r3[] = {0.0, -2.12159572323E5, 2.30661510616E5,
		      2.74647644705E4, -4.02621119975E4,
		      -2.29660729780E3, -1.16328495004E5,
		      -1.46025937511E5, -2.42357409629E4,
		      -5.70691009324E2};

static const double r4[] = {0.0, 2.79195317918525E-1, 4.917317610505968E-1,
		      6.92910599291889E-2, 3.350343815022304E0,
		      6.012459259764103E0};

static const double alr2pi = 9.18938533204673E-1;
static const double xlge = 5.10E6;

/***********************************************************************
 *
 * FUNCTION:     aln_gamma
 * 
 * DESCRIPTION:  Compute ln(gamma(x)) using the AS 245 algortithm
 *
 * PARAMETERS:   x
 *
 * RETURN:       ln(gamma(x))
 *      
 ***********************************************************************/
static double aln_gamma(double x) SPECFUN;
static double aln_gamma(double x)
{
	double alngam,y,x1,x2;

	if (x >= 1E305)
		return NaN;
	if (x <= 0)
		return NaN;
	
	alngam = 0.0;

	/* Calculation for 0 < X < 0.5 and 0.5 <= X < 1.5 combined */
	if (x < 1.5) {
		if (x < 0.5) {
			alngam = -log(x);
			y = x + 1.0;
			
			/* Test whether X < machine epsilon */
			if (y == 1.0)
				return alngam;
		} else {
			alngam = 0.0;
			y = x;
			x = (x - 0.5) - 0.5;
		}
		alngam = alngam + x * ((((r1[5]*y + r1[4])*y + r1[3])*y
				+ r1[2])*y + r1[1]) / ((((y+r1[9])*y 
				 + r1[8])*y + r1[7])*y+r1[6]);
		return alngam;
	}
	/* Calculation for 1.5 <= X < 4.0 */
	if (x < 4.0) {
		y = (x - 1.0) - 1.0;
		alngam = y * ((((r2[5]*x + r2[4])*x + r2[3])*x + r2[2])*x
			      + r2[1]) / ((((x + r2[9])*x + r2[8])*x 
					   + r2[7])*x + r2[6]);
		return alngam;
	}
	/* Calculation for 4.0 <= X < 12.0 */
	if (x < 12.0) {
		alngam = ((((r3[5]*x + r3[4])*x + r3[3])*x + r3[2])*x + r3[1]) /
			((((x + r3[9])*x + r3[8])*x + r3[7])*x + r3[6]);
		return alngam;
	}
	/* Calculation for X >= 12.0 */
	y = log(x);
	alngam = x * (y - 1.0) - y/2.0 + alr2pi;
	if (x > xlge)
		return alngam;
	x1 = 1.0 / x;
	x2 = x1 * x1;
	alngam = alngam + x1 * ((r4[3]*x2 + r4[2])*x2 + r4[1]) /
		((x2 + r4[5])*x2 + r4[4]);
	return alngam;
}

#define MAX_ITER 200

/***********************************************************************
 *
 * FUNCTION:     gser
 * 
 * DESCRIPTION:  Compute incomplete gamma function, power series
 *
 * PARAMETERS:   a, x
 *
 * RETURN:       P(a,x)
 *      
 ***********************************************************************/
static double gser(double a, double x) SPECFUN;
static double gser(double a, double x)
{
  double gln, ap, sum, del;
  int n;

  gln = aln_gamma(a);
  if (x == 0)
    return(0);
  ap = a;
  sum = 1.0/a;
  del = sum;
  for (n=0; n<MAX_ITER; n++) {
    ap++;
    del*=x/ap;
    sum+=del;
    if (fabs(del < fabs(sum)*EPSILON))
      break;
  }
  /* Screw it if the series failed to converge after the maximum number
     of iterations.  We need to determine a way of signalling to the user
     that the iteration failed to converge--next version maybe. */
  return(sum*exp(-x + a*log(x) - gln));
}

/***********************************************************************
 *
 * FUNCTION:     gcf
 * 
 * DESCRIPTION:  Compute complementary incomplete gamma function,
 *               continued fraction rep.
 *
 * PARAMETERS:   a, x
 *
 * RETURN:       Q(a,x)
 *      
 ***********************************************************************/
static double gcf(double a, double x) SPECFUN;
static double gcf(double a, double x)
{
  double gln, g_old, g, a0, a1, b0, b1, fac, ana, anf;
  int n;

  gln = aln_gamma(a);
  g_old = 0;
  a0 = 1;
  a1 = x;
  b0 = 0;
  b1 = 1;
  fac = 1;
  g = 1;
  for (n=1; n<=MAX_ITER; n++) {
    ana = (double)n - a;
    a0 = (a1 + a0*ana)*fac;
    b0 = (b1 + b0*ana)*fac;
    anf = n*fac;
    a1 = x*a0+anf*a1;
    b1 = x*b0+anf*b1;
    if (a1 != 0.0) {
      fac = 1.0/a1;
      g = b1*fac;
      if (fabs((g-g_old)/g) < EPSILON)
	break;
      g_old = g;
    }
  }
  return(exp(-x + a*log(x) - gln)*g);
}

/***********************************************************************
 *
 * FUNCTION:     incgamma
 * 
 * DESCRIPTION:  Compute incomplete gamma function
 *
 * PARAMETERS:   a, x
 *
 * RETURN:       P(a,x)
 *      
 ***********************************************************************/
static double incgamma(double a, double x) SPECFUN;
static double incgamma(double a, double x)
{
  if (x < a+1)
    return(gser(a, x));
  return(1.0-gcf(a, x));
}

/***********************************************************************
 *
 * FUNCTION:     erf
 * 
 * DESCRIPTION:  Compute error function
 *
 * PARAMETERS:   x
 *
 * RETURN:       erf(x)
 *      
 ***********************************************************************/
static double erf(double x) SPECFUN;
static double erf(double x)
{
  if (x < 0)
    return(-incgamma(0.5, x*x));
  return(incgamma(0.5, x*x));
}

/***********************************************************************
 *
 * FUNCTION:     erf
 * 
 * DESCRIPTION:  Compute complementary error function
 *
 * PARAMETERS:   x
 *
 * RETURN:       erfc(x)
 *      
 ***********************************************************************/
static double erfc(double x) SPECFUN;
static double erfc(double x)
{
	double xx;
	
	xx = x*x;
	if (x < 0) {
		if (xx < 1.5)
			return 1.0 + gser(0.5, xx);
		else
			return 2.0 - gcf(0.5, xx);
	} else {
		if (xx < 1.5)
			return 1.0 - gser(0.5, xx);
		else
			return gcf(0.5, xx);
	}
}

/***********************************************************************
 *
 * FUNCTION:     betacf
 * 
 * DESCRIPTION:  Compute incomplete beta function, continued fraction
 *
 * PARAMETERS:   a, b, x
 *
 * RETURN:       ibeta(a, b, x)
 *      
 ***********************************************************************/
static double betacf(double a, double b, double x) SPECFUN;
static double betacf(double a, double b, double x)
{
  double am, bm, az, qab, qap, qam, bz, em, tem, ap, bp, d, app, bpp, aold;
  int m;

  am = bm = az = 1;
  qab = a+b;
  qap = a+1;
  qam = a-1;
  bz = 1-qab*x/qap;
  for (m=1; m<=MAX_ITER; m++) {
    em = m;
    tem = em + em;
    d = em*(b-m)*x/((qam+tem)*(a+tem));
    ap = az + d*am;
    bp = bz + d*bm;
    d = -(a+em)*(qab+em)*x/((a+tem)*(qap+tem));
    app = ap + d*az;
    bpp = bp + d*bz;
    aold = az;
    am = ap/bpp;
    bm = bp/bpp;
    az = app/bpp;
    bz = 1.0;
    if (fabs(az - aold) < EPSILON*fabs(az))
      break;
  }
  return(az);
}

/***********************************************************************
 *
 * FUNCTION:     ibeta
 * 
 * DESCRIPTION:  Compute incomplete beta function
 *
 * PARAMETERS:   a, b, x
 *
 * RETURN:       ibeta(a, b, x)
 *      
 ***********************************************************************/
static double incbeta(double a, double b, double x) SPECFUN;
static double incbeta(double a, double b, double x)
{
  double bt;

  bt = (x == 0 || x == 1) ? 0 :
    exp(aln_gamma(a+b) - aln_gamma(a) -
	    aln_gamma(b) + a*log(x) + b*log(1-x));
  if (x < (a+1)/(a+b+2))		/* use continued fraction directly */
    return(bt*betacf(a, b, x)/a);
  /* use continued fraction after making symmetry transform */
  return(1-bt*betacf(b, a, 1-x)/b);
}

/***********************************************************************
 *
 * FUNCTION:     math_gambeta
 * 
 * DESCRIPTION:  Compute Gamma, Beta, and related functions
 *
 * PARAMETERS:   On stack - up to three numbers, depending on function
 *
 * RETURN:       On stack - 1 number
 *      
 ***********************************************************************/
CError
math_gambeta(Functype *func, CodeStack *stack)
{
  double z, w, v, res;
  CError err;

  err = stack_get_val(stack, &z, real);
  if (err)
    return(err);

  switch (func->num) {
  case MATH_BETA:
    err=stack_get_val(stack, &w, real);
    if (err)
      return(err);
    if (z < 0 || w < 0)
      return(c_badarg);
    v = z+w;
    res = exp(aln_gamma(z) + aln_gamma(w) - aln_gamma(z+w));
    break;
  case MATH_IGAMMA:
    err=stack_get_val(stack, &w, real);
    if (err)
      return(err);
    if (z < 0 || w <= 0)
      return(c_badarg);
    res = incgamma(w, z);
    break;
  case MATH_ERF:
    res = erf(z);
    break;
  case MATH_ERFC:
    res = erfc(z);
    break;
  case MATH_IBETA:
    err = stack_get_val2(stack, &w, &v, real);
    if (err)
      return(err);
    if (z < 0 || z > 1 || w <= 0 || v <= 0)
      return(c_badarg);
    res = incbeta(w, v, z);
    break;
  default:
    return(c_internal);
  }
  return stack_add_val(stack, &res, real);
}

/***********************************************************************
 *
 * FUNCTION:     polevl
 * 
 * DESCRIPTION:  Evaluate a polynomial given an array of its coefficients

 * PARAMETERS:   x, coeffs, degree N
 *
 * RETURN:       value of polynomial at x
 *      
 ***********************************************************************/
double polevl(double x, double *coef, int N )
{
  double p;
  int i;

  p = coef[N-1];
  for (i=N-2; i>=0; i--)
    p = p*x + coef[i];
  return(p);
}


Index: history.h
===================================================================
RCS file: /cvsroot/easycalc/PPCport/core/mlib/history.h,v
retrieving revision 1.2
retrieving revision 1.3
diff -C2 -d -r1.2 -r1.3
*** history.h	24 Jun 2009 21:39:59 -0000	1.2
--- history.h	17 Oct 2009 13:48:34 -0000	1.3
***************
*** 33,37 ****
  	union {
  		dbStackItem item;
! 		TCHAR *text;
  	}u;
  }tHistory;
--- 33,37 ----
  	union {
  		dbStackItem item;
! 		TCHAR text[1];
  	}u;
  }tHistory;

--- NEW FILE: matrix.h ---
/*
 *   $Id: matrix.h,v 1.1 2009/10/17 13:48:34 mapibid Exp $
 *
 *   Scientific Calculator for Palms.
 *   Copyright (C) 1999,2000,2001 Ondrej Palkovsky
 *
 *   This program is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation; either version 2 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program; if not, write to the Free Software
 *   Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *   Boston, MA  02110-1301, USA.
 * 
 *  You can contact me at 'on...@pe...'.
*/

#ifndef _MATRIX_H_
#define _MATRIX_H_

#define MATRIX_MAX  50

#define MATRIX_INPUT   2001
#define MATRIX_DET     2002
#define MATRIX_ECHELON 2003
#define MATRIX_IDENTITY 2004
#define MATRIX_QRQ 2005
#define MATRIX_QRR 2006
#define MATRIX_QRS 2007

CError matrix_input(Functype *func,CodeStack *stack) MATFUNC;
void matrix_delete(Matrix *m) MATFUNC;
Matrix * matrix_new(Int16 rows,Int16 cols) MATFUNC;
Matrix * matrix_dup(Matrix *m) MATFUNC;
CError matrix_oper(Functype *func,CodeStack *stack) MATFUNC;
CError matrix_oper_real(Functype *func,CodeStack *stack) MATFUNC;
CError matrix_func(Functype *func,CodeStack *stack) MATFUNC;
CError matrix_func2(Functype *func,CodeStack *stack) MATFUNC;
CError matrix_item(Functype *func, CodeStack *stack) MATFUNC;
CError  matrix_dim(Functype *func, CodeStack *stack) MATFUNC;
CError matrix_func3(Functype *func,CodeStack *stack) MATFUNC;

#endif

--- NEW FILE: elliptic.cpp ---
/*
 *   $Id: elliptic.cpp,v 1.1 2009/10/17 13:48:34 mapibid Exp $
 * 
 *   Scientific Calculator for Palms.
 *   Copyright (C) 1999,2000 Ondrej Palkovsky
 *   This file Copyright (C) 2000 Rafael R. Sevilla
 *
 *   This program is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation; either version 2 of the License, or
 *   (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program; if not, write to the Free Software
 *   Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *   Boston, MA  02110-1301, USA.
 * 
 *  You can contact Ondrej at 'on...@pe...'.
 *  You can contact me at 'di...@pa...'.
 *
 *  Portions of this code are from the Cephes Math Library, Release 2.8
 *  Copyright (C) 1984, 1987, 1989, 2000 by Stephen L. Moshier
 */
#include "StdAfx.h"

#include "compat/PalmOS.h"
#include "compat/MathLib.h"
#include "konvert.h"
#include "funcs.h"
#include "mathem.h"
#include "stack.h"
#include "core/prefs.h"
#include "display.h"
#include "specfun.h"

#define PIO2 (M_PIl*0.5)
#define C1 1.3862943611198906188E0

/***********************************************************************
 *
 * FUNCTION:     ellc1
 * 
 * DESCRIPTION:  Compute complete elliptic integral of the first kind
 *
 * PARAMETERS:   eccentricity m
 *
 * RETURN:       K(m)
 *      
 ***********************************************************************/
static double ellc1(double x) SPECFUN;
static double ellc1(double x)
{
  static double P[] =
  {
    1.38629436111989062502E0,
    9.65735902811690126535E-2,
    3.08851465246711995998E-2,
    1.49380448916805252718E-2,
    8.79078273952743772254E-3,
    6.18901033637687613229E-3,
    6.87489687449949877925E-3,
    9.85821379021226008714E-3,
    7.97404013220415179367E-3,
    2.28025724005875567385E-3,
    1.37982864606273237150E-4
  };

  static double Q[] =
  {
    4.99999999999999999821E-1,
    1.24999999999870820058E-1,
    7.03124996963957469739E-2,
    4.88280347570998239232E-2,
    3.73774314173823228969E-2,
    3.01204715227604046988E-2,
    2.39089602715924892727E-2,
    1.54850516649762399335E-2,
    5.94058303753167793257E-3,
    2.94078955048598507511E-5,
  };

  if (x > EPSILON)
    return(polevl(x, P, 10) - log(x)*polevl(x, Q, 10));
  return(C1 - 0.5*log(x));
}

/***********************************************************************
 *
 * FUNCTION:     elli1
 * 
 * DESCRIPTION:  Compute incomplete elliptic integral of the first kind
 *
 * PARAMETERS:   eccentricity m and modulus phi
 *
 * RETURN:       F(m, phi)
 *      
 ***********************************************************************/
static double elli1(double phi, double m) SPECFUN;
static double elli1(double phi, double m)
{
  double a, b, c, e, temp, t, K;
  int d, mod, sign, npio2;

  if( m == 0.0 )
    return( phi );
  a = 1.0 - m;
  if( a == 0.0 ) {
    return(log(tan((PIO2 + phi)/2.0)));
  }
  npio2 = (int) floor(phi/PIO2);
  if(npio2 & 1)
    npio2 += 1;
  if(npio2) {
    K = ellc1(a);
    phi = phi - npio2 * PIO2;
  } else
    K = 0.0;

  if( phi < 0.0 ) {
    phi = -phi;
    sign = -1;
  } else
    sign = 0;
  b = sqrt(a);
  t = tan( phi );
  if( fabs(t) > 10.0 ) {
    /* Transform the amplitude */
    e = 1.0/(b*t);
    /* ... but avoid multiple recursions.  */
    if( fabs(e) < 10.0 ) {
      e = atan(e);
      if(npio2 == 0)
	K = ellc1(a);
      temp = K - elli1(e, m);
      goto done;
    }
  }
  a = 1.0;
  c = sqrt(m);
  d = 1;
  mod = 0;

  while(fabs(c/a) > EPSILON) {
    temp = b/a;
    phi = phi + atan(t*temp) + mod * M_PIl;
    mod = (int) ((phi + PIO2)/M_PIl);
    t = t * ( 1.0 + temp )/( 1.0 - temp * t * t );
    c = ( a - b )/2.0;
    temp = sqrt( a * b );
    a = ( a + b )/2.0;
    b = temp;
    d += d;
  }

  temp = (atan(t) + mod * M_PIl)/(d * a);

 done:
  if( sign < 0 )
    temp = -temp;
  temp += npio2 * K;
  return(temp);
}

/***********************************************************************
 *
 * FUNCTION:     ellc2
 * 
 * DESCRIPTION:  Compute complete elliptic integral of the second kind
 *
 * PARAMETERS:   eccentricity m
 *
 * RETURN:       E(m)
 *      
 ***********************************************************************/
static double ellc2(double m) SPECFUN;
static double ellc2(double m)
{
  static double P[] = {
    1.00000000000000000299E0,
    4.43147180560990850618E-1,
    5.68051945617860553470E-2,
    2.18317996015557253103E-2,
    1.15688436810574127319E-2,
    7.58395289413514708519E-3,
    7.77395492516787092951E-3,
    1.07350949056076193403E-2,
    8.68786816565889628429E-3,
    2.50888492163602060990E-3,
    1.53552577301013293365E-4,
  };

  static double Q[] = {
    2.49999999999888314361E-1,
    9.37499997197644278445E-2,
    5.85936634471101055642E-2,
    4.27180926518931511717E-2,
    3.34833904888224918614E-2,
    2.61769742454493659583E-2,
    1.68862163993311317300E-2,
    6.50609489976927491433E-3,
    1.00962792679356715133E-3,
    3.27954898576485872656E-5,
  };

  if (m == 0)
    return(1.0);
  return(polevl(m, P, 10) - log(m)*(m*polevl(m, Q, 9)));
}

/***********************************************************************
 *
 * FUNCTION:     elli2
 * 
 * DESCRIPTION:  Compute incomplete elliptic integral of the second kind
 *
 * PARAMETERS:   eccentricity m and modulus phi
 *
 * RETURN:       F(m, phi)
 *      
 ***********************************************************************/
static double elli2(double phi, double m) SPECFUN;
static double elli2(double phi, double m)
{
  double a, b, c, e, temp;
  double lphi, t, E;
  int d, mod, npio2, sign;

  if( m == 0.0 )
    return( phi );
  lphi = phi;
  npio2 = (int) floor( lphi/PIO2 );
  if( npio2 & 1 )
    npio2 += 1;
  lphi = lphi - npio2 * PIO2;
  if( lphi < 0.0 ) {
    lphi = -lphi;
    sign = -1;
  } else {
    sign = 1;
  }
  a = 1.0 - m;
  E = ellc2(a);
  if( a == 0.0 ) {
    temp = sin(lphi);
    goto done;
  }
  t = tan(lphi);
  b = sqrt(a);
  if( fabs(t) > 10.0 ) {
    /* Transform the amplitude */
    e = 1.0/(b*t);
    /* ... but avoid multiple recursions.  */
    if( fabs(e) < 10.0 ) {
      e = atan(e);
      temp = E + m * sin( lphi ) * sin( e ) - elli2( e, m );
      goto done;
    }
  }
  c = sqrt(m);
  a = 1.0;
  d = 1;
  e = 0.0;
  mod = 0;

  while( fabs(c/a) > EPSILON ) {
    temp = b/a;
    lphi = lphi + atan(t*temp) + mod * M_PIl;
    mod = (int) ((lphi + PIO2)/M_PIl);
    t = t * ( 1.0 + temp )/( 1.0 - temp * t * t );
    c = ( a - b )/...
 
[truncated message content]