700 lines
16 KiB
C
700 lines
16 KiB
C
// Symbolic expressions are built by connecting U structs.
|
|
//
|
|
// For example, (a b + c) is built like this:
|
|
//
|
|
// _______ _______ _______
|
|
// |CONS |--->|CONS |----------------------------->|CONS |
|
|
// | | | | | |
|
|
// |_______| |_______| |_______|
|
|
// | | |
|
|
// ___v___ ___v___ _______ _______ ___v___
|
|
// |ADD | |CONS |--->|CONS |--->|CONS | |SYM c |
|
|
// | | | | | | | | | |
|
|
// |_______| |_______| |_______| |_______| |_______|
|
|
// | | |
|
|
// ___v___ ___v___ ___v___
|
|
// |MUL | |SYM a | |SYM b |
|
|
// | | | | | |
|
|
// |_______| |_______| |_______|
|
|
|
|
typedef struct U {
|
|
union {
|
|
struct {
|
|
struct U *car; // pointing down
|
|
struct U *cdr; // pointing right
|
|
} cons;
|
|
struct {
|
|
struct U *binding; // symbol's value binding
|
|
struct U *binding2; // symbol's function binding
|
|
} sym;
|
|
char *str;
|
|
struct tensor *tensor;
|
|
struct {
|
|
unsigned int *a, *b; // rational number a over b
|
|
} q;
|
|
double d;
|
|
} u;
|
|
unsigned char k, tag;
|
|
} U;
|
|
|
|
// the following enum is for struct U, member k
|
|
|
|
enum {
|
|
CONS,
|
|
NUM,
|
|
DOUBLE,
|
|
STR,
|
|
TENSOR,
|
|
SYM,
|
|
};
|
|
|
|
// the following enum is for indexing the symbol table
|
|
|
|
enum {
|
|
|
|
// standard functions first, then nil, then everything else
|
|
|
|
ABS,
|
|
ADD,
|
|
ADJ,
|
|
AND,
|
|
ARCCOS,
|
|
ARCCOSH,
|
|
ARCSIN,
|
|
ARCSINH,
|
|
ARCTAN,
|
|
ARCTANH,
|
|
ARG,
|
|
ATOMIZE,
|
|
BESSELJ,
|
|
BESSELY,
|
|
BINDING2,
|
|
BINOMIAL,
|
|
BREAK,
|
|
CARAC,
|
|
CEILING,
|
|
CHARPOLY,
|
|
CHECK,
|
|
CLEAR,
|
|
CLS,
|
|
COEFF,
|
|
CONDENSE,
|
|
CONJ,
|
|
CONTRACT,
|
|
CONVOLUTION,
|
|
COS,
|
|
COSH,
|
|
DEGREE,
|
|
DENOMINATOR,
|
|
DERIVATIVE,
|
|
DET,
|
|
DIM,
|
|
DIRAC,
|
|
DISPLAY,
|
|
DIVISORS,
|
|
DO,
|
|
DOT,
|
|
DRAW,
|
|
DSOLVE,
|
|
EIGEN,
|
|
EIGENVAL,
|
|
EIGENVEC,
|
|
ERF,
|
|
ERFC,
|
|
EVAL,
|
|
EXP,
|
|
EXPAND,
|
|
EXPCOS,
|
|
EXPSIN,
|
|
FACTOR,
|
|
FACTORIAL,
|
|
FACTORPOLY,
|
|
FILTER,
|
|
FLOATF,
|
|
FLOOR,
|
|
FOR,
|
|
FOURIER,
|
|
GAMMA,
|
|
GCD,
|
|
HEAVISIDE,
|
|
HERMITE,
|
|
HILBERT,
|
|
IMAG,
|
|
INDEX,
|
|
INNER,
|
|
INTEGRAL,
|
|
INV,
|
|
INVFOURIER,
|
|
INVG,
|
|
ISINTEGER,
|
|
ISPRIME,
|
|
LAGUERRE,
|
|
LCM,
|
|
LEGENDRE,
|
|
LOG,
|
|
MAG,
|
|
MOD,
|
|
MULTIPLY,
|
|
NOT,
|
|
NUMERATOR,
|
|
OPERATOR,
|
|
OR,
|
|
OUTER,
|
|
POWER,
|
|
PRIME,
|
|
PRINT,
|
|
PRODUCT,
|
|
PROG,
|
|
QUOTE,
|
|
RANK,
|
|
RATIONALIZE,
|
|
REAL,
|
|
YYRECT,
|
|
RETURN,
|
|
ROOTS,
|
|
SETQ,
|
|
SGN,
|
|
SIMFAC,
|
|
SIMPLIFY,
|
|
SIN,
|
|
SINH,
|
|
SQRT,
|
|
STOP,
|
|
SUBST,
|
|
SUM,
|
|
TAB,
|
|
TAN,
|
|
TANH,
|
|
TAYLOR,
|
|
TCHEBYCHEVT,
|
|
TCHEBYCHEVU,
|
|
TEST,
|
|
TESTEQ,
|
|
TESTGE,
|
|
TESTGT,
|
|
TESTLE,
|
|
TESTLT,
|
|
TRACE,
|
|
TRANSPOSE,
|
|
UNIT,
|
|
WEDGE,
|
|
ZERO,
|
|
|
|
NIL, // nil goes here, after standard functions
|
|
|
|
AUTOEXPAND,
|
|
EXPOMODE,
|
|
LAST,
|
|
TTY,
|
|
YYVOID,
|
|
YYE,
|
|
YYLAST,
|
|
YYTRUE,
|
|
YYFALSE,
|
|
|
|
// symbols appearing above are printed in roman type
|
|
|
|
// the following symbols are printed in italic type
|
|
|
|
PI,
|
|
SYMBOL_A,
|
|
SYMBOL_B,
|
|
SYMBOL_C,
|
|
SYMBOL_D,
|
|
SYMBOL_N,
|
|
SYMBOL_R,
|
|
SYMBOL_T,
|
|
SYMBOL_X,
|
|
SYMBOL_Y,
|
|
SYMBOL_Z,
|
|
|
|
USR_SYMBOLS, // this must be last
|
|
};
|
|
|
|
#define E YYE
|
|
|
|
#define TOS 1000000
|
|
#define BUF 10000
|
|
|
|
#define MAX_PROGRAM_SIZE 100001
|
|
#define MAXPRIMETAB 10000
|
|
|
|
#define _USE_MATH_DEFINES // for MS C++
|
|
|
|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <ctype.h>
|
|
#include <fcntl.h>
|
|
#include <string.h>
|
|
#include <setjmp.h>
|
|
#include <math.h>
|
|
#include <errno.h>
|
|
|
|
#define MAXDIM 24
|
|
|
|
typedef struct tensor {
|
|
int ndim;
|
|
int dim[MAXDIM];
|
|
int nelem;
|
|
U *elem[1];
|
|
} T;
|
|
|
|
struct display {
|
|
int h, w, n;
|
|
struct {
|
|
int c, x, y;
|
|
} a[1]; // a for array
|
|
};
|
|
|
|
struct text_metric {
|
|
int ascent, descent, width;
|
|
};
|
|
|
|
extern U **frame;
|
|
|
|
#define iscons(p) ((p)->k == CONS)
|
|
#define isrational(p) ((p)->k == NUM)
|
|
#define isdouble(p) ((p)->k == DOUBLE)
|
|
#define isnum(p) (isrational(p) || isdouble(p))
|
|
#define isstr(p) ((p)->k == STR)
|
|
#define istensor(p) ((p)->k == TENSOR)
|
|
#define issymbol(p) ((p)->k == SYM)
|
|
#define iskeyword(p) (issymbol(p) && symbol_index(p) < NIL)
|
|
|
|
#define car(p) (iscons(p) ? (p)->u.cons.car : nil)
|
|
#define cdr(p) (iscons(p) ? (p)->u.cons.cdr : nil)
|
|
#define caar(p) car(car(p))
|
|
#define cadr(p) car(cdr(p))
|
|
#define cdar(p) cdr(car(p))
|
|
#define cddr(p) cdr(cdr(p))
|
|
#define caadr(p) car(car(cdr(p)))
|
|
#define caddr(p) car(cdr(cdr(p)))
|
|
#define cadar(p) car(cdr(car(p)))
|
|
#define cdadr(p) cdr(car(cdr(p)))
|
|
#define cddar(p) cdr(cdr(car(p)))
|
|
#define cdddr(p) cdr(cdr(cdr(p)))
|
|
#define caaddr(p) car(car(cdr(cdr(p))))
|
|
#define cadadr(p) car(cdr(car(cdr(p))))
|
|
#define caddar(p) car(cdr(cdr(car(p))))
|
|
#define cdaddr(p) cdr(car(cdr(cdr(p))))
|
|
#define cadddr(p) car(cdr(cdr(cdr(p))))
|
|
#define cddddr(p) cdr(cdr(cdr(cdr(p))))
|
|
#define caddddr(p) car(cdr(cdr(cdr(cdr(p)))))
|
|
#define cadaddr(p) car(cdr(car(cdr(cdr(p)))))
|
|
#define cddaddr(p) cdr(cdr(car(cdr(cdr(p)))))
|
|
#define caddadr(p) car(cdr(cdr(car(cdr(p)))))
|
|
#define cdddaddr(p) cdr(cdr(cdr(car(cdr(cdr(p))))))
|
|
#define caddaddr(p) car(cdr(cdr(car(cdr(cdr(p))))))
|
|
|
|
//#define arg1(p) car(cdr(p))
|
|
//#define arg2(p) car(cdr(cdr(p)))
|
|
//#define arg3(p) car(cdr(cdr(cdr(p))))
|
|
|
|
#define isadd(p) (car(p) == symbol(ADD))
|
|
#define ispower(p) (car(p) == symbol(POWER))
|
|
#define isfactorial(p) (car(p) == symbol(FACTORIAL))
|
|
|
|
//#define numerator mp_numerator
|
|
void mp_numerator(void);
|
|
|
|
//#define denominator mp_denominator
|
|
void mp_denominator(void);
|
|
|
|
#define isfraction mp_isfraction
|
|
int mp_isfraction(U *);
|
|
|
|
#define scan_integer bignum_scan_integer
|
|
void bignum_scan_integer(char *);
|
|
|
|
#define scan_float bignum_scan_float
|
|
void bignum_scan_float(char *);
|
|
|
|
void bignum_power_number(int);
|
|
|
|
void bignum_truncate(void);
|
|
|
|
void mp_set_bit(unsigned int *, unsigned int);
|
|
void mp_clr_bit(unsigned int *, unsigned int);
|
|
|
|
extern U *pop();
|
|
|
|
extern int tos;
|
|
extern int expanding;
|
|
extern int expomode;
|
|
extern int conjugating;
|
|
extern int fmt_x;
|
|
extern int fmt_index;
|
|
extern int fmt_level;
|
|
extern int verbosing;
|
|
extern int floating;
|
|
extern char program_buf[];
|
|
extern int primetab[MAXPRIMETAB];
|
|
extern int esc_flag;
|
|
extern U *nil;
|
|
extern U *stack[];
|
|
|
|
extern U *p1, *p2, *p3, *p4, *p5, *p6, *p7, *p8;
|
|
extern U *formal_arg[6];
|
|
extern U *tmp;
|
|
extern U *nil;
|
|
extern U *zero, *one, *imaginaryunit;
|
|
extern U *table_of_integrals;
|
|
extern U *table_of_fourier;
|
|
extern U *meta_a;
|
|
extern U *meta_b;
|
|
extern U *meta_c;
|
|
extern U *meta_n;
|
|
extern U *meta_x;
|
|
|
|
extern U *alloc();
|
|
extern U *alloc_tensor();
|
|
extern U *inner_product(U *, U *);
|
|
extern U *number();
|
|
extern U *pop();
|
|
|
|
extern int cmp(U **, U**);
|
|
|
|
extern U *copy_number(U *);
|
|
double convert_rational_to_double(U *);
|
|
|
|
extern void run(char *);
|
|
extern void save(void);
|
|
extern void restore(void);
|
|
extern void d_tensor_tensor(void);
|
|
extern void d_tensor_scalar(void);
|
|
extern void d_scalar_tensor(void);
|
|
extern void d_scalar_scalar(void);
|
|
extern int equal(U *, U *);
|
|
extern int equaln(U *, int);
|
|
extern int length(U *);
|
|
extern void push(U *);
|
|
extern void addk(int);
|
|
extern void multiply_all(int);
|
|
extern void divide(void);
|
|
extern void multiply(void);
|
|
extern void logarithm(void);
|
|
extern void add(void);
|
|
extern int lessp(U *, U *);
|
|
extern void list(int);
|
|
extern void sine(void);
|
|
extern void cosine(void);
|
|
extern void expcos(void);
|
|
extern void expsin(void);
|
|
extern void negate(void);
|
|
extern void negate_expand(void);
|
|
extern void negate_noexpand(void);
|
|
extern void init(void);
|
|
extern void defn(void);
|
|
extern void test(char *, char **, int);
|
|
extern void test_all(void);
|
|
extern void read_file(char *);
|
|
extern int scan(char *);
|
|
extern void print(U *);
|
|
extern void qfree(U *);
|
|
extern void init(void);
|
|
extern U *number(int);
|
|
extern void push_double(double);
|
|
extern void push_rational(int, int);
|
|
extern void out_of_memory(void);
|
|
extern void maxdim_error(void);
|
|
extern void cons(void);
|
|
extern int isfraction(U *);
|
|
extern void subtract(void);
|
|
extern int ispoly(U *, U *);
|
|
extern int ispoly_expr(U *, U *);
|
|
extern void subst(void);
|
|
extern int coeff(void);
|
|
extern void eval(void);
|
|
extern void eval_predicate(void);
|
|
extern void eval_noexpand(void);
|
|
extern void evalat(void);
|
|
extern int find(U *, U *);
|
|
extern int ispoly_term(U *, U *);
|
|
extern int ispoly_factor(U *, U *);
|
|
extern int isposint(U *);
|
|
extern int compare_tensors(U *, U *);
|
|
extern void inv(void);
|
|
extern void inverse(void);
|
|
extern void term_plus_term(void);
|
|
extern void print_number(U *);
|
|
extern void tensor_plus_tensor(void);
|
|
extern void factor_times_factor(int *);
|
|
extern void push_number(int);
|
|
extern void push_integer(int);
|
|
extern int pop_integer(void);
|
|
extern void tensor_times_scalar(void);
|
|
extern void scalar_times_tensor(void);
|
|
extern void spow(void);
|
|
extern int is_numeric_factor(U *);
|
|
extern void multiply_rational_powers(int);
|
|
extern void qnumer(void);
|
|
extern void qdenom(void);
|
|
extern int qfit(U *);
|
|
extern void qinv(void);
|
|
extern void qtrunc(void);
|
|
extern void qpow(void);
|
|
extern void conj_tensor(void);
|
|
extern void build_tensor(int);
|
|
extern U *alloc_tensor(int);
|
|
extern void print_unsigned_term(U *);
|
|
extern void print_unsigned_factor(U *);
|
|
extern void print_factor(U *);
|
|
extern void print_tensor(U *);
|
|
extern void print_tensor_inner(U *, int, int *);
|
|
extern void print_str(char *);
|
|
extern void print_char(int);
|
|
extern void print_num(U *, int);
|
|
extern void print_double(U *, int flag);
|
|
extern void print_function_definition(U *);
|
|
extern void print_arg_list(U *);
|
|
extern void print_lisp(U *);
|
|
extern void print1(U *);
|
|
extern void eval(void);
|
|
extern void iterate(int);
|
|
extern void adj(void);
|
|
extern void arccos(void);
|
|
extern void arcsin(void);
|
|
extern void arctan(void);
|
|
extern void contract(void);
|
|
extern void scot(void);
|
|
extern void scsc(void);
|
|
extern void cross(void);
|
|
extern void curl(void);
|
|
extern void derivative(void);
|
|
extern int integer(void);
|
|
extern void det(void);
|
|
extern void inner(void);
|
|
extern void inv(void);
|
|
extern void outer(void);
|
|
extern void ssec(void);
|
|
extern void tangent(void);
|
|
extern void taylor(void);
|
|
extern void trace(void);
|
|
extern void transpose(void);
|
|
extern void conjugate(void);
|
|
extern void tensor_plus_tensor(void);
|
|
extern void determinant(int);
|
|
extern int is_square_matrix(U *);
|
|
extern void cofactor(U *, int, int, int);
|
|
extern void component(int);
|
|
extern void push_cars(U *);
|
|
extern void mp_test(void);
|
|
extern void push_int(int);
|
|
extern void integral(void);
|
|
extern int cmp_expr(U *, U *);
|
|
extern void dsolve(void);
|
|
extern void sexp(void);
|
|
extern void swap(void);
|
|
extern void distill(void);
|
|
extern void distill_sum(void);
|
|
extern void distill_product(void);
|
|
extern void scanlisp(char *);
|
|
extern int distilly();
|
|
extern void simplify(void);
|
|
extern void simfac(void);
|
|
extern void peek(void);
|
|
extern void peek2(void);
|
|
extern void expand(void);
|
|
extern void expand1(void);
|
|
extern void mul2(void);
|
|
extern void factorial(void);
|
|
extern void hilbert(void);
|
|
extern void swap(void);
|
|
extern void guess(void);
|
|
extern int isnegativeterm(U *);
|
|
char *get_name(U *);
|
|
int new_name(char *);
|
|
void init_alloc(void);
|
|
void detg(void);
|
|
void invg(void);
|
|
void eval_tensor(void);
|
|
void test_a(void);
|
|
void test_b(void);
|
|
void test_c(void);
|
|
void test_d(void);
|
|
void test_h(void);
|
|
void test_i(void);
|
|
void test_k(void);
|
|
void testp(void);
|
|
void testp_prime(void);
|
|
void testq(void);
|
|
void testq_prime(void);
|
|
void gc(void);
|
|
void push_frame(int);
|
|
void pop_frame(int);
|
|
void binomial(void);
|
|
void charpoly(void);
|
|
void init_stack(void);
|
|
void push_zero_matrix(int, int);
|
|
void push_identity_matrix(int);
|
|
void push_symbol(int);
|
|
char *get_printname(U *);
|
|
U *symbol(int);
|
|
void power(void);
|
|
void exponential(void);
|
|
void invert_matrix(void);
|
|
void power_tensor(void);
|
|
void init_integers(void);
|
|
void stop(char *);
|
|
void top_eval(void);
|
|
void hermite(void);
|
|
void prime(void);
|
|
int iszero(U *);
|
|
int isplusone(U *);
|
|
int isminusone(U *);
|
|
int isinteger(U *);
|
|
int isintegerfactor(U *);
|
|
int isnonnegativeinteger(U *);
|
|
int iseveninteger(U *);
|
|
void ssinh(void);
|
|
void scosh(void);
|
|
void gcd(void);
|
|
void laguerre(void);
|
|
void lcm(void);
|
|
void legendre(void);
|
|
void test_binomial(void);
|
|
void test_hermite(void);
|
|
void test_laguerre(void);
|
|
void test_legendre(void);
|
|
void square(void);
|
|
void rationalize(void);
|
|
int isnegativenumber(U *);
|
|
int isnegative(U *);
|
|
void test_gcd(void);
|
|
void factor(void);
|
|
void test_factor(void);
|
|
void fmt_char(int);
|
|
void fmt_fraction(int, int, int);
|
|
void fmt_test(U *);
|
|
void scan2(char *);
|
|
void test_is_prime(void);
|
|
void test_quickfactor(void);
|
|
void factor_number(void);
|
|
void factor_poly(void);
|
|
void degree(void);
|
|
void divisors(void);
|
|
void divisors_onstack(void);
|
|
void factor_small_number(void);
|
|
int sign(int);
|
|
void factorpoly(void);
|
|
void multiply_all_noexpand(int);
|
|
void vectorize(int);
|
|
void quickfactor(void);
|
|
void alloc_qnum(void);
|
|
void restore_frame(U **);
|
|
void multiply_noexpand(void);
|
|
void absval(void);
|
|
void variables(void);
|
|
void coeff_cooked(void);
|
|
void display(U *);
|
|
void print_out(void);
|
|
void prog(void);
|
|
void prog_return(void);
|
|
void loop(void);
|
|
void loop_break(void);
|
|
void __print(void);
|
|
double pop_double(void);
|
|
void bignum_float(void);
|
|
void new_string(char *);
|
|
void printline(U *);
|
|
void tokenize(char *);
|
|
int scan_next(void);
|
|
void wedge2(void);
|
|
void wedge3(void);
|
|
int save_symbols(int);
|
|
void restore_symbols(int);
|
|
void printstr(char *);
|
|
void printchar(int);
|
|
void printchar_nowrap(int);
|
|
void condense(void);
|
|
extern void ysinh(void);
|
|
extern void ycosh(void);
|
|
extern char logbuf[];
|
|
extern void logout(char *);
|
|
extern void errout(void);
|
|
extern void clear(void);
|
|
extern void clear_symbols(void);
|
|
extern void clear_term(void);
|
|
extern void numerator(void);
|
|
extern void denominator(void);
|
|
extern void reciprocate(void);
|
|
|
|
#define MSIGN(p) (((int *) (p))[-2])
|
|
#define MLENGTH(p) (((int *) (p))[-1])
|
|
|
|
#define MZERO(p) (MLENGTH(p) == 1 && (p)[0] == 0)
|
|
#define MEQUAL(p, n) (MLENGTH(p) == 1 && (long long) MSIGN(p) * (p)[0] == (n))
|
|
|
|
extern unsigned int *mnew(int);
|
|
extern void mfree(unsigned int *);
|
|
extern unsigned int *mint(int);
|
|
extern unsigned int *mcopy(unsigned int *);
|
|
extern int mcmp(unsigned int *, unsigned int *);
|
|
extern int mcmpint(unsigned int *, int);
|
|
extern unsigned int *madd(unsigned int *, unsigned int *);
|
|
extern unsigned int *msub(unsigned int *, unsigned int *);
|
|
extern unsigned int *mmul(unsigned int *, unsigned int *);
|
|
extern unsigned int *mdiv(unsigned int *, unsigned int *);
|
|
extern unsigned int *mmod(unsigned int *, unsigned int *);
|
|
void mdivrem(unsigned int **, unsigned int **, unsigned int *, unsigned int *);
|
|
extern unsigned int *mpow(unsigned int *, unsigned int);
|
|
extern unsigned int *mroot(unsigned int *, unsigned int);
|
|
extern unsigned int *msqrt(unsigned int *);
|
|
extern unsigned int *mgcd(unsigned int *, unsigned int *);
|
|
extern unsigned int *mmodpow(unsigned int *, unsigned int *, unsigned int *);
|
|
extern void mshiftright(unsigned int *);
|
|
extern unsigned int *mscan(char *);
|
|
extern char *mstr(unsigned int *);
|
|
extern int mprime(unsigned int *);
|
|
extern int mtotal;
|
|
extern int issymbolic(U *);
|
|
|
|
extern void qadd(void);
|
|
extern void qsub(void);
|
|
extern void qmul(void);
|
|
extern void qdiv(void);
|
|
extern void qpow(void);
|
|
|
|
// these functions handle mixed-mode arithmetic (rationals and doubles)
|
|
|
|
extern int compare_numbers(U *, U *);
|
|
extern void add_numbers(void);
|
|
extern void subtract_numbers(void);
|
|
extern void multiply_numbers(void);
|
|
extern void divide_numbers(void);
|
|
extern void power_numbers(void);
|
|
extern void negate_number(void);
|
|
extern void invert_number(void);
|
|
extern void gcd_numbers(void);
|
|
|
|
extern void besselj(void);
|
|
extern void bessely(void);
|
|
extern void carac(void);
|
|
extern void convolution(void);
|
|
extern void dirac(void);
|
|
extern void gamma(void);
|
|
extern void yerf(void);
|
|
extern void yerfc(void);
|
|
extern void fourier(void);
|
|
extern void heaviside(void);
|
|
extern void invfourier(void);
|
|
extern void sgn(void);
|
|
extern void tchebychevT(void);
|
|
extern void tchebychevU(void);
|
|
extern double erfc(double);
|
|
extern void std_symbol(char *, int);
|
|
extern U *usr_symbol(char *);
|
|
extern int symbol_index(U *);
|
|
extern int iscomplexnumber(U *);
|
|
extern void real(void);
|
|
extern void imag(void);
|
|
extern void mag(void);
|
|
extern void arg(void);
|
|
extern void rect(void);
|
|
extern void simplify_trig(void);
|
|
extern int equalq(U *, int, int);
|
|
|
|
// In case I forget to use YY...
|
|
|
|
#undef TRUE
|
|
#undef FALSE
|