erf stuff

This commit is contained in:
George Weigt 2005-08-05 12:28:02 -07:00
parent f717e70734
commit cb4cb6463f
12 changed files with 108 additions and 132 deletions

View file

@ -538,41 +538,13 @@ static struct {
{"Paste\tCtrl+V", ID_PASTE},
{0, 0},
{"*", 0},
{"*", 0},
{"About", ID_ABOUT},
{"Memory", ID_MEMORY},
{"Copy display to clipboard", ID_COPY_DISPLAY},
{"Create script from command history", ID_CREATE_SCRIPT},
{0, 0},
#if 0
{"Help", 0},
{"Help Pages", ID_HELP_PAGES},
{"-", 0},
{"Type ^ for exponent", ID_HELP_EXPONENT},
{"Type a space to multiply", ID_HELP_MULTIPLY},
{"How to draw a graph", ID_HELP_DRAW},
{"How to factor a polynomial", ID_HELP_FACTOR_POLYNOMIAL},
{"How to factor a number", ID_HELP_FACTOR_NUMBER},
{"How to define a symbol", ID_HELP_SYMBOL},
{"How to define a function", ID_HELP_FUNCTION},
//{"A special note about functions", ID_HELP_SPECIAL_NOTE},
{"How to define a vector", ID_HELP_TYPE_VECTOR},
{"How to define a matrix", ID_HELP_TYPE_MATRIX},
{"How to multiply a matrix and vector", ID_HELP_MATRIX_TIMES_VECTOR},
{"How to invert a matrix", ID_HELP_INVERT_MATRIX},
{"How to draw a parametric graph", ID_HELP_DRAW_CIRCLE},
{"Sample Scripts", 0},
{"Gamma Matrix Algebra", ID_SAMPLE_GMA},
{"Vector Calculus", ID_SAMPLE_VC},
{"Rotation Matrix", ID_SAMPLE_RM},
{"Quantum Harmonic Oscillator", ID_SAMPLE_QHO},
{"Hydrogenic Wavefunctions", ID_SAMPLE_HW},
{"Static Spherical Metric", ID_SAMPLE_SSM},
{"Free Particle Dirac Equation", ID_SAMPLE_FPDE},
{0, 0},
//{"About", ID_HELP_ABOUT},
{0, 0},
#endif
{0, 0},
};
@ -723,6 +695,7 @@ LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
// case ID_PRINT:
// do_print();
// break;
case ID_UNDO:
if (running)
break;
@ -755,11 +728,14 @@ LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
else
SendMessage(edit_window, WM_PASTE, 0, 0);
break;
// "asterisk" pull-down menu
case ID_ABOUT:
if (running)
break;
goto_calc_mode();
printstr("Version 106 Help is available at eigenmath.com\n");
printstr("Version 106 - Help is available at eigenmath.com\n");
update_display();
break;
case ID_MEMORY:
@ -777,58 +753,9 @@ LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
break;
do_create_script();
break;
#if 0
case ID_HELP_EXPONENT:
do_main_help(1);
break;
case ID_HELP_MULTIPLY:
do_main_help(2);
break;
case ID_HELP_DRAW:
do_main_help(3);
break;
case ID_HELP_FACTOR_POLYNOMIAL:
do_main_help(4);
break;
case ID_HELP_FACTOR_NUMBER:
do_main_help(5);
break;
case ID_HELP_SYMBOL:
do_main_help(6);
break;
case ID_HELP_FUNCTION:
do_main_help(7);
break;
case ID_HELP_SPECIAL_NOTE:
do_main_help(8);
break;
case ID_HELP_TYPE_VECTOR:
do_main_help(9);
break;
case ID_HELP_TYPE_MATRIX:
do_main_help(10);
break;
case ID_HELP_MATRIX_TIMES_VECTOR:
do_main_help(11);
break;
case ID_HELP_INVERT_MATRIX:
do_main_help(12);
break;
case ID_HELP_DRAW_CIRCLE:
do_main_help(13);
break;
case ID_HELP_PAGES:
HtmlHelp(main_window, "help.chm", HH_DISPLAY_TOPIC, NULL);
break;
#endif
case ID_EDIT_SCRIPT:
if (running)
break;
if (edit_mode == 0)
goto_edit_mode();
else
goto_calc_mode();
break;
// window buttons, clear, draw, etc.
case ID_CLEAR:
do_button("clear");
break;
@ -847,32 +774,18 @@ LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
case ID_INTEGRAL:
do_button("integral");
break;
case ID_EDIT_SCRIPT:
if (running)
break;
if (edit_mode == 0)
goto_edit_mode();
else
goto_calc_mode();
break;
case ID_RUN_SCRIPT:
run_script();
break;
#if 0
case ID_SAMPLE_GMA:
do_example(0);
break;
case ID_SAMPLE_VC:
do_example(1);
break;
case ID_SAMPLE_RM:
do_example(2);
break;
case ID_SAMPLE_QHO:
do_example(3);
break;
case ID_SAMPLE_HW:
do_example(4);
break;
case ID_SAMPLE_SSM:
do_example(5);
break;
case ID_SAMPLE_FPDE:
do_example(6);
break;
#endif
default:
return DefWindowProc(hWnd, message, wParam, lParam);
}

View file

@ -88,8 +88,6 @@ gc(void)
untag(p7);
untag(p8);
untag(last);
for (i = 0; i < tos; i++)
untag(stack[i]);

View file

@ -9,7 +9,7 @@ U *p5;
U *p6;
U *p7;
U *p8;
U *last;
U *tmp;
U *nil;
U *formal_arg[6];

7
defs.h
View file

@ -205,6 +205,10 @@ enum {
SYMBOL_Y,
SYMBOL_Z,
TTY,
// system symbols
YYLAST,
};
#define MAXDIM 24
@ -309,7 +313,6 @@ extern U *stack[];
extern U *p1, *p2, *p3, *p4, *p5, *p6, *p7, *p8;
extern U *formal_arg[6];
extern U *tmp;
extern U *last;
extern U *nil;
extern U *zero, *one, *imaginaryunit;
extern U *table_of_integrals;
@ -655,4 +658,4 @@ extern void invfourier(void);
extern void sgn(void);
extern void tchebychevT(void);
extern void tchebychevU(void);
extern double erfc(double);

View file

@ -898,6 +898,11 @@ static char *s[] = {
"d(g(f(x)),x)",
"d(g(f(x)),x)",
// other functions
"d(erf(x))-2*exp(-x^2)/sqrt(pi)",
"0",
};
void

View file

@ -55,8 +55,18 @@ static int draw_count;
void
eval_draw(void)
{
p2 = caddr(p1);
p1 = cadr(p1);
// save "last" so we can do, for example,
//
// x^2
// draw
// integral
//
// and have integral do integral(x^2)
push(symbol(YYLAST)->u.sym.binding);
p2 = caddr(p1); // free variable in function
p1 = cadr(p1); // function to draw
// if first arg is a symbol then evaluate it
@ -79,12 +89,15 @@ eval_draw(void)
p2 = symbol(SYMBOL_X);
}
push(p1);
push(p2);
push(p1); // function to draw
push(p2); // free variable in function
draw();
push(nil);
symbol(YYLAST)->u.sym.binding = pop();
push(nil); // so no result is printed
// also, "last" is not modified when result is "nil"
}
static void

17
erf.cpp
View file

@ -30,16 +30,16 @@ yerf(void)
static void
yyerf(void)
{
//double d;
double d;
p1 = pop();
#if 0
if (p1->k == DOUBLE) {
d = erf(p1->u.d);
d = 1.0 - erfc(p1->u.d);
push_double(d);
return;
}
#endif
if (isnegativeterm(p1)) {
push_symbol(ERF);
push(p1);
@ -59,10 +59,13 @@ static char *s[] = {
"erf(a)",
"erf(a)",
"erf(0.0) + 1", // add 1 to round off
"1",
"float(erf(0)) + 1", // add 1 to round off
"1",
#if 0
"float(erf(0))",
"0",
"float(erf(1))",
"0.842701",
#endif

View file

@ -3,6 +3,8 @@
// Author : philippe.billet@noos.fr
//
// erfc(x)
//
// GW Added erfc() from Numerical Recipes in C
//
//-----------------------------------------------------------------------------
@ -29,30 +31,51 @@ yerfc(void)
static void
yyerfc(void)
{
//double d;
double d;
p1 = pop();
#if 0
if (p1->k == DOUBLE) {
d = erfc(p1->u.d);
push_double(d);
return;
}
#endif
push_symbol(ERFC);
push(p1);
list(2);
return;
}
// from Numerical Recipes in C
#ifndef LINUX
double
erfc(double x)
{
double t, z, ans;
z = fabs(x);
t = 1.0 / (1.0 + 0.5 * z);
ans=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
t*(-0.82215223+t*0.17087277)))))))));
return x >= 0.0 ? ans : 2.0-ans;
}
#endif
static char *s[] = {
"erfc(a)",
"erfc(a)",
#if 0
"erfc(0.0)",
"1",
"float(erfc(0))",
"1",
#if 0
"float(erfc(1))",
"0.157299",
#endif

View file

@ -974,7 +974,7 @@ eval(void)
// save last non-nil result
if (stack[tos - 1] != nil)
last = stack[tos - 1];
symbol(YYLAST)->u.sym.binding = stack[tos - 1];
}
static void

View file

@ -24,7 +24,6 @@ init(void)
p7 = nil;
p8 = nil;
last = nil;
varlist = nil;
table_of_integrals = nil;
@ -153,7 +152,7 @@ init(void)
define_symbol("wedge", WEDGE);
define_symbol("zero", ZERO);
// built-in symbols
// user symbols
define_variable("autoexpand", AUTOEXPAND);
define_variable("~exp", E); // tilde so sort puts it after scalar symbols
@ -174,6 +173,10 @@ init(void)
define_variable("y", SYMBOL_Y);
define_variable("z", SYMBOL_Z);
// system symbols
define_variable("$last", YYLAST);
meta_a = get_symbol("$a");
meta_b = get_symbol("$b");
meta_c = get_symbol("$c");

View file

@ -934,6 +934,10 @@ char *integrals[] = {
"exp(a*x^2)*exp(b)*(x^2/a-1/(a^2))/2",
NULL,
"exp(a*x^2)",
"-i*sqrt(pi)*erf(i*sqrt(a)*x)/sqrt(a)/2",
NULL,
NULL,
};
@ -1581,6 +1585,17 @@ static char *s[] = {
"#572",
"integral(cosh(X)^2,X)-sinh(2*X)/4-X/2",
"0",
// test integral(exp(a*x^2))
"integral(exp(a*x^2))+i*sqrt(pi)*erf(i*sqrt(a)*x)/sqrt(a)/2",
"0",
"integral(exp(-x^2))-sqrt(pi)*erf(x)/2",
"0",
"integral(exp(-3*x^2))-sqrt(pi/3)*erf(sqrt(3)*x)/2",
"0",
};
void

View file

@ -73,7 +73,7 @@ run(char *s)
p2 = pop();
check_stack();
symbol(LAST)->u.sym.binding = last;
symbol(LAST)->u.sym.binding = symbol(YYLAST)->u.sym.binding;
symbol(LAST)->u.sym.binding2 = nil;
// print string w/o quotes