eigenmath/fourier.cpp

1100 lines
19 KiB
C++
Raw Permalink Normal View History

2005-07-30 21:37:29 +02:00
//-----------------------------------------------------------------------------
//
// Author : philippe.billet@noos.fr
//
// Fourier transform :
// The formalism used is : integral(f(t)*exp(-i*t*x)) dt
// The following properties are implemented :
// - linearity :fourier(a*f(t)+b*g(t),t)=a*fourier(f(t),t)+b*fourier(g(t),t)
// - derivation : fourier(d(f(t),t),t)=i*t*fourier(f(t),t)
// - product by t : fourier(t*f(t),t)=-i*d(fourier(f(t),t),t)
2006-01-06 21:26:36 +01:00
// - integral : fourier(integral(f(t),t),t)=fourier(f(t),t)/(i*t)
2005-07-30 21:37:29 +02:00
// - division by t : fourier(f(t)/t,t)=integral(fourier(f(t),t),t)/i
// - product of convolution : fourier(convolution(f(t),g(t)),t)=fourier(f(t),t)*fourier(g(t),t)
// - fourier of fourier : fourier(fourier(f(t),t),t)=2*pi*f(-t)
//
// The formalism used is : integral(f(t)*exp(-i*t*x)) dt
// N.B.: the is no change of variable : fourier(x,x)=2*i*pi*d(dirac(x),x)
//-----------------------------------------------------------------------------
#include "stdafx.h"
#include "defs.h"
static void yfourier(void);
static void fourier_of_sum(void);
static void fourier_of_product(void);
static void fourier_of_convolution(void);
static void fourier_of_derivative1(void);
static void fourier_of_derivative2(void);
static void fourier_of_integral1(void);
static void fourier_of_integral2(void);
static void fourier_of_fourier(void);
static void fourier_of_invfourier(void);
static void fourier_of_form(void);
2006-01-06 21:26:36 +01:00
static void fourier_of_inverse_trig(void);
2005-07-30 21:37:29 +02:00
static int match(U *, U *, U *, int, int);
static void scan_fouriers(void);
2006-09-20 18:09:25 +02:00
static int exp_flag;
2006-01-06 21:26:36 +01:00
2005-07-30 21:37:29 +02:00
void
eval_fourier(void)
{
2006-01-06 21:26:36 +01:00
int expotemp;
2006-05-06 01:27:26 +02:00
expotemp=exp_flag;
2006-01-06 21:26:36 +01:00
2006-05-06 01:27:26 +02:00
exp_flag=0;
2006-01-06 21:26:36 +01:00
push(cadr(p1));
eval();
2006-01-16 20:37:31 +01:00
if (caddr(p1) == symbol(NIL))
2006-01-06 21:26:36 +01:00
guess();
else
push(caddr(p1));
fourier();
p5=pop();
// print(p5);
// printf("\n");
// printf("%d\n",length(p5));
2006-05-06 01:27:26 +02:00
exp_flag=1;
2005-07-30 21:37:29 +02:00
push(cadr(p1));
eval();
2006-01-16 20:37:31 +01:00
if (caddr(p1) == symbol(NIL))
2005-07-30 21:37:29 +02:00
guess();
2005-10-04 02:02:16 +02:00
else
2005-07-30 21:37:29 +02:00
push(caddr(p1));
fourier();
2006-01-06 21:26:36 +01:00
p6=pop();
// print(p6);
// printf("\n");
// printf("%d\n",length(p6));
//
if (find(p5,symbol(FOURIER))==1)
push(p6);
else
if ((length(p6) <= length(p5)) || (find(p6,symbol(FOURIER))==1))
push(p5);
else
push(p6);
2006-05-06 01:27:26 +02:00
exp_flag=expotemp;
2005-07-30 21:37:29 +02:00
}
void
fourier(void)
{
save();
yfourier();
restore();
}
static void
yfourier(void)
{
2006-01-06 21:26:36 +01:00
2006-01-16 20:37:31 +01:00
if (table_of_fourier == symbol(NIL))
2005-07-30 21:37:29 +02:00
scan_fouriers();
p2 = pop();
p1 = pop();
push(caddr(p1));
negate();
p3 = pop();
2006-01-06 21:26:36 +01:00
2005-07-30 21:37:29 +02:00
if (car(p1) == symbol(ADD))
fourier_of_sum();
else if (car(p1) == symbol(MULTIPLY))
fourier_of_product();
else if (car(p1) == symbol(CONVOLUTION))
fourier_of_convolution();
else if (car(p1) == symbol(DERIVATIVE)){
if (caddr(p1) == p2)
fourier_of_derivative1();
else
fourier_of_derivative2();}
else if (car(p1) == symbol(INTEGRAL)){
if (caddr(p1) == p2)
fourier_of_integral1();
else
fourier_of_integral2();}
else if (car(p1) == symbol(FOURIER) && (caddr(p1) == p2 || p3 == p2))
fourier_of_fourier();
else if (car(p1) == symbol(INVFOURIER) && caddr(p1) == p2)
fourier_of_invfourier();
2006-01-06 21:26:36 +01:00
else if ( find(p1,symbol(ARCCOS))==1 ||
find(p1,symbol(ARCCOSH))==1 ||
find(p1,symbol(ARCSIN))==1 ||
find(p1,symbol(ARCSINH))==1 ||
find(p1,symbol(ARCTAN))==1 ||
find(p1,symbol(ARCTANH))==1/* ||
find(p1,symbol(HEAVISIDE))==1 ||
find(p1,symbol(SGN))==1*/)
fourier_of_inverse_trig();
2005-07-30 21:37:29 +02:00
else
fourier_of_form();
}
2006-01-06 21:26:36 +01:00
2005-07-30 21:37:29 +02:00
static void
fourier_of_form(void)
{
int h;
// save meta symbols
push(meta_a->u.sym.binding);
push(meta_b->u.sym.binding);
push(meta_c->u.sym.binding);
push(meta_n->u.sym.binding);
push(meta_x->u.sym.binding);
meta_x->u.sym.binding = p2;
h = tos;
push(one);
push(p1);
push(p2);
distill();
p3 = table_of_fourier;
while (iscons(p3)) {
p4 = car(p3);
if (match(p1, car(p4), cddr(p4), h, tos))
break;
p3 = cdr(p3);
}
if (iscons(p3)) {
push(cadr(p4));
eval();
} else {
push_symbol(FOURIER);
push(p1);
push(p2);
list(3);
}
p3 = pop();
tos = h;
// restore meta symbols
meta_x->u.sym.binding = pop();
meta_n->u.sym.binding = pop();
meta_c->u.sym.binding = pop();
meta_b->u.sym.binding = pop();
meta_a->u.sym.binding = pop();
push(p3);
}
static void
fourier_of_sum(void)
{
p1 = cdr(p1);
push(car(p1));
push(p2);
fourier();
p1 = cdr(p1);
while (iscons(p1)) {
push(car(p1));
push(p2);
fourier();
add();
p1 = cdr(p1);
}
}
static void
fourier_of_derivative1(void)
{
p1 = cdr(p1);
push(imaginaryunit);
push(p2);
multiply();
push(car(p1));
push(p2);
fourier();
multiply();
}
static void
fourier_of_derivative2(void)
{
push(cadr(p1));
push(p2);
fourier();
push(caddr(p1));
derivative();
}
static void
fourier_of_integral1(void)
{
p1 = cdr(p1);
push(car(p1));
push(p2);
fourier();
2006-01-06 21:26:36 +01:00
p3=pop();
/* print(p3);
printstr("\n");*/
push(p3);
2005-07-30 21:37:29 +02:00
push(imaginaryunit);
push(p2);
multiply();
divide();
2006-01-06 21:26:36 +01:00
// if ((find(p3,symbol(FOURIER))==0) &&
// (find(p3,symbol(INTEGRAL))==0) &&
// (find(p3,symbol(DERIVATIVE))==0)) {
// push_symbol(PI);
// push(p3);
// push(p2);
// push_integer(0);
// evalat();
// multiply();
// add();
// }
2005-07-30 21:37:29 +02:00
}
static void
fourier_of_integral2(void)
{
push(cadr(p1));
push(p2);
fourier();
push(caddr(p1));
integral();
}
static void
fourier_of_fourier(void)
{
p1 = cdr(p1);
push_symbol(PI);
push_integer(2);
multiply();
push(car(p1));
push(p2);
push(p2);
negate();
subst();
multiply();
}
static void
fourier_of_invfourier(void)
{
p1 = cdr(p1);
push(car(p1));
}
static void
fourier_of_convolution(void)
{ p1 = cdr(p1);
push(car(p1));
push(p2);
fourier();
p1 = cdr(p1);
while (iscons(p1)) {
push(car(p1));
push(p2);
fourier();
multiply();
p1 = cdr(p1);
}
2006-01-06 21:26:36 +01:00
}
2005-07-30 21:37:29 +02:00
2006-01-06 21:26:36 +01:00
static void
fourier_of_inverse_trig(void)
{
push(p1);
push(p2);
derivative();
push(p2);
fourier();
push(imaginaryunit);
push(p2);
multiply();
divide();
2005-07-30 21:37:29 +02:00
}
static void
fourier_of_product(void)
{
2006-01-06 21:26:36 +01:00
int h1, h2,n=0;
2005-07-30 21:37:29 +02:00
h1 = tos;
p3 = cdr(p1);
while (iscons(p3)) {
if (!find(car(p3), p2))
push(car(p3));
p3 = cdr(p3);
}
h2 = tos;
p3 = cdr(p1);
2006-01-06 21:26:36 +01:00
2005-07-30 21:37:29 +02:00
while (iscons(p3)) {
2006-01-06 21:26:36 +01:00
if (find(car(p3), p2)){
push(car(p3));}
if (caar(p3) == symbol(POWER) && cadar(p3) == p2 && isinteger(caddar(p3))){
push(caddar(p3));
n = pop_integer();
}
if (car(p3) == p2){
push_integer(1);
n = pop_integer();
}
2005-07-30 21:37:29 +02:00
p3 = cdr(p3);
}
if (tos - h2 == 0){
push_integer(2);
push_symbol(PI);
multiply();
push(p2);
dirac();
2006-01-06 21:26:36 +01:00
multiply();}
2005-07-30 21:37:29 +02:00
else {
multiply_all(tos - h2);
p1 = pop();
2006-01-06 21:26:36 +01:00
if (!(n==0)){
push(p1);
push(p1);
p3=pop();
2005-07-30 21:37:29 +02:00
push(p2);
2006-01-06 21:26:36 +01:00
push_integer(n);
power();
divide();
p1=pop();
}
if (!isplusone(p1)) {
fourier_of_form();
if (!(n==0)) {
if (n > 0) {
while(n > 0){
push(p2);
derivative();
push(imaginaryunit);
multiply();
n=n-1;}
2005-07-30 21:37:29 +02:00
}
2006-01-06 21:26:36 +01:00
else {
n=-n;
while(n > 0){
push(p2);
integral();
push(imaginaryunit);
push_integer(-1);
multiply();
multiply();
n=n-1;}
2005-07-30 21:37:29 +02:00
}
2006-01-06 21:26:36 +01:00
}
}
else {
push(p3);
push(p2);
fourier();
}
2005-07-30 21:37:29 +02:00
}
multiply_all(tos - h1);
}
static int
match(U *actual, U *formal, U *caveats, int h1, int h2)
{
int i, j, k, l;
save();
for (i = h1; i < h2; i++) {
for (j = h1; j < h2; j++) {
for (k = h1; k < h2; k++) {
for (l = h1; l < h2; l++) {
meta_a->u.sym.binding = stack[i];
meta_b->u.sym.binding = stack[j];
meta_c->u.sym.binding = stack[k];
meta_n->u.sym.binding = stack[l];
// check caveats
p1 = caveats;
while (iscons(p1)) {
push(car(p1));
eval();
p2 = pop();
2006-01-11 00:30:19 +01:00
if (iszero(p2))
2005-07-30 21:37:29 +02:00
break;
p1 = cdr(p1);
}
if (iscons(p1))
continue;
// actual == formal?
push(formal);
eval();
if (equal(actual, pop())) {
restore();
return 1;
}
}
}
}
}
restore();
return 0;
}
2006-01-06 21:26:36 +01:00
char * fouriers[] = {
2005-07-30 21:37:29 +02:00
"a",
"2*pi*a*dirac(x)",
NULL,
2006-01-06 21:26:36 +01:00
"dirac(a*x)",
"a",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"dirac(a*x+b)",
"exp(i*b*x/a)/abs(a)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"heaviside(x)",
"pi*dirac(x)-i/x",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"heaviside(x+a)",
"pi*dirac(x)-i*exp(i*a*x)/x",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"heaviside(x)*exp(a*x)",
"(-a-i*x)/(a^2+x^2)",
"isnegativeterm(a)==1",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"sgn(a*x)",
"-2*i*sgn(a)/x",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"sgn(a*x+b)",
"-2*i*exp(i*b*x/a)*sgn(a)/x",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*x)*sgn(x)",
"-2*i/(x+i*a)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"carac(x,a,b)",
"i*(exp(-i*b*x)-exp(-i*a*x))/x",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"carac(x,-a,a)*exp(b*x)",
"2*sin(a*(x+i*b))/(x+i*b)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"carac(x,-a,a)*abs(x)",
"2*a*(sin(a*x)/(a*x)-2*(sin(a*x/2)/(a*x))^2)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"sgn(x)*carac(x,-a,a)",
"-4*i*sin(a*x/2)*sin(a*x/2)/x",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"shah(a*x)",
"2*pi*shah(2*pi*x/a)/a",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"shah(a*x+b)",
"2*pi*exp(i*b*x/a)*shah(2*pi*x/a)/a",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"besselj(a*x,0)",
"2/(abs(a)*sqrt(1-(x/a)^2))",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"besselj(a*x+b,0)",
"2*exp(i*b*x/a)/(abs(a)*sqrt(1-(x/a)^2))",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*x)*besselj(b*x,0)",
"2/(abs(b)*sqrt(1-((x+i*a)/b)^2))",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"1/sqrt(a+b*x^2)",
"pi*besselj(-sgn(b/a)*x*sqrt(-a/b),0)/sqrt(a*abs(b))",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"1/sqrt(1+a*(x+b)^2)",
"pi*besselj(-sgn(a)x*sqrt(-1/a),0)*exp(i*b*x)/sqrt(abs(a))",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*x)/sqrt(1+b*x^2)",
"pi*besselj(-sgn(b)(-i*x+i*a)*sqrt(-1/b),0)/sqrt(abs(b))",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"besselj(a*x,1)/(2*a*x)",
"sqrt(1-(x/a)^2)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"sqrt(1+a*x^2)",
"pi*besselj(i*sqrt(a)*x,1)/(x*abs(i*sqrt(a)))",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"sqrt(1+a*(x+b)^2)",
"pi*besselj(i*sqrt(a)*x,1)*exp(i*b*x)/(x*abs(i*sqrt(a)))",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*x)*sqrt(1+b*x^2)",
"pi*besselj(i*sqrt(b)*(x+i*a),1)/((x+i*a)*abs(i*sqrt(b)))",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"erf(a*x)",
"-i*2*exp(-(x/(2*a))^2)/x",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"erf(a*x+b)",
"-i*2*exp(-(x/(2*a))^2)*exp(i*b*x/a)/x",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"x^a",
"2*(i^a)*pi*d(dirac(x),x,a)",
"and(a>0,isinteger(a)==1)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"(x+a)^b",
"2*(i^b)*pi*d(dirac(x),x,b))*exp(i*a*x)",
"and(b>0,isinteger(b)==1)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"x^a",
"(-1)^(-a)*i^(-a)*x^(-a-1)*pi*sgn(x)/(-a-1)!",
"and(a<0,isinteger(a)==1)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"(x+a)^b",
"(-1)^(-b)*i^(-b)*x^(-b-1)*exp(i*a*x)*pi*sgn(x)/(-b-1)!",
"and(b<0,isinteger(b)==1)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
2005-07-30 21:37:29 +02:00
"1/(a*x+b)",
"-i*exp(i*b*x/a)*sgn(x)*pi/a",
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*x)/(c*x+b)",
"-i*exp(i*b*(x+i*a)/c)*sgn(x+i*a)*pi/c",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"1/(a*x^2+b)",
"exp(-sqrt(b/a)*abs(x))*pi/(sqrt(b*a))",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"1/(a*x+b)^2",
"-abs(x/a)*exp(i*b*x/a)*pi/abs(a)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*x)/(b*x^2+c)",
"exp(-sqrt(c/b)*abs(x+i*a))*pi/(sqrt(b*c))",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"1/(a*(x+b)^2+c)",
"exp(i*b*x)*exp(sqrt(c/a)*abs(x))*pi/(sqrt(c*a))",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"1/(a*x^2+b*x)",
"(-i*pi*sgn(x)+i*pi*exp(i*b*x/a)*sgn(x))/b",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"1/abs(a*x)",
"(-2*euler-2*log(abs(x/a)))/abs(a)",
NULL,
"1/abs(a*x+b)",
"(-2*euler-2*log(abs(x/a)))*exp(i*b*x/a)/abs(a)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*x)/abs(x)",
"-2*euler-2*log(abs(x+i*a))",
NULL,
"exp(a*x)/abs(b*x+c)",
"(-2*euler-2*log(abs((x+i*a)/b)))*exp(i*c*(x+i*a)/b)/abs(b)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"log(a*abs(x)))",
"(-2*pi*euler*dirac(x/abs(a))-pi/abs(x/a))/abs(a)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"log(abs(a*x+b))",
"(-2*pi*euler*dirac(x/abs(a))-pi/abs(x/a))*exp(i*b*x/a)/abs(a)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*x)*log(abs(b*x))",
"(-2*pi*euler*dirac((x+i*a)/b)-pi/abs((x+i*a)/b))/abs(b)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*x)*log(abs(b*x+c))",
"(-2*pi*euler*dirac((x+i*a)/b)-pi/abs(x/b))*exp(i*c*(x+i*a)/b)/abs(b)",
2005-07-30 21:37:29 +02:00
NULL,
"abs(a*x)",
"-2/(abs(a)*(x/a)^2)",
NULL,
"abs(a*x+b)",
"-2*exp(i*b*x/a)/(abs(a)*(x/a)^2)",
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*x)*abs(b*x)",
"-2/(abs(b)*((x+i*a)/b)^2)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"abs(x)^a",
"-2*Gamma(a+1)*sin(pi*a/2)*abs(x)^(-a-1)",
"and(not(abs(a)==1),isinteger(a)==0)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*x)*abs(x)^b",
"-2*Gamma(b+1)*sin(pi*b/2)*abs((x+i*a))^(-b-1)",
"and(not(abs(b)==1),isinteger(b)==0)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"abs(x+a)^b",
"-2*exp(i*a*x)*Gamma(b+1)*sin(pi*b/2)*abs(x)^(-b-1)",
"and(not(abs(b)==1),isinteger(b)==0)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*x)*abs(x+c)^b",
"-2*exp(i*c*(x+i*a))*Gamma(b+1)*sin(pi*b/2)*abs((x+i*a))^(-b-1)",
"and(not(abs(b)==1),isinteger(b)==0)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"abs(x)^a*sgn(x)",
"-2*i*Gamma(a+1)*cos(pi*a/2)*sgn(x)*abs(x)^(-a-1)",
"and(not(abs(a)==1),isinteger(a)==0)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*x)*abs(x)^b*sgn(x)",
"-2*i*Gamma(b+1)*cos(pi*b/2)*sgn((x+i*a))*abs(x)^(-b-1)",
"and(not(abs(b)==1),isinteger(b)==0)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"abs(x+b)^a*sgn(x+b)",
"-2*i*exp(i*b*x)*Gamma(a+1)*cos(pi*a/2)*sgn(x)*abs(x)^(-a-1)",
"and(not(abs(a)==1),isinteger(a)==0)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"abs(x)^a*heaviside(x)",
"-i*abs(x)^(-a-1)*Gamma(1+a)*(cos(a*pi/2)*sgn(x)-i*sin(a*pi/2))",
"and(not(abs(a)==1),isinteger(a)==0)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*x)*abs(x)^b*heaviside(x)",
"-i*abs((x+i*a))^(-b-1)*Gamma(1+b)*(cos(b*pi/2)*sgn((x+i*a))-i*sin(b*pi/2))",
"and(not(abs(b)==1),isinteger(b)==0)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"abs(x+a)^b*heaviside(x+a)",
"-i*exp(i*a*x)*abs(x)^(-b-1)*Gamma(1+b)*(cos(b*pi/2)*sgn(x)-i*sin(b*pi/2))",
"and(not(abs(b)==1),isinteger(b)==0)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*x)*abs(x+c)^b*heaviside(x+c)",
"-i*exp(i*c*(x+i*a))*abs((x+i*a))^(-b-1)*Gamma(1+b)*(cos(b*pi/2)*sgn((x+i*a))-i*sin(b*pi/2))",
"and(not(abs(b)==1),isinteger(b)==0)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*x)",
"2*pi*dirac(x+i*a)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*x+b*x)",
"2*pi*dirac(x+i*a+i*b)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*x+b)",
"2*pi*dirac(x+i*a)*exp(b)",
NULL,
"heaviside(x)*exp(a*x)",
"(-a-i*x)/(a^2+x^2)",
"isnegativeterm(a)==1",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*x)*heaviside(x)",
"(i*a+i*x)/(a^2+x^2)+pi*(dirac(x+i*a)-dirac(x-i*a))",
"and(real(a)==0,not(imag(a)==0))",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"heaviside(x+b)*exp(a*(x+b))",
"(a+x)*exp(i*b*x)/(a^2+x^2)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"heaviside(x)*exp(a*x+b*x)",
"1/(-b+(i*a+x)*i)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*abs(x))",
"-2*a/(a^2+x^2)",
"isnegativeterm(a)==1",
NULL,
"exp(a*abs(b*x+c))",
"-2*a*exp(i*c*x/b)/((a^2+(x/b)^2)*abs(b))",
"isnegativeterm(a)==1",
NULL,
"exp(a*x+b*abs(x))",
"-2*b/(b^2+(x+i*a)^2)",
NULL,
"2*a/(a^2+(x+b)^2)",
"-exp(-i*b*x-a*abs(x))",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"sgn(x)*exp(a*abs(x))",
"-2*i*x/(a^2+x^2)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"sgn(x+b)*exp(a*abs(b*x+c))",
"-2*i*(x/b)*exp(i*c*x/b)/(a^2+(x/b)^2)",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*x^2)",
"sqrt(pi/(2*abs(imag(a))))*exp(-i*x^2/(4*imag(a)))*(1+i*sgn(imag(a)))",
"and(real(a)==0,not(imag(a)==0))",
2005-07-30 21:37:29 +02:00
NULL,
2006-01-06 21:26:36 +01:00
"exp(a*(x+b)^2)",
"sqrt(pi/(2*abs(imag(a))))*exp(-i*x^2/(4*imag(a)imag(a)))*(1+i*sgn(imag(a)))*exp(i*b*x)",
"and(real(a)==0,not(imag(a)==0))",
NULL,
"exp(a*x^2)",
"sqrt(pi/abs(a))*exp(x^2/(4*a))",
"and(imag(a)==0,isnegativeterm(a)==1)",
NULL,
"exp(a*(x+b)^2)",
"sqrt(pi/abs(a))*exp(x^2/(4*a))*exp(i*b*x)",
"and(imag(a)==0,isnegativeterm(a)==1)",
2005-07-30 21:37:29 +02:00
NULL,
};
static void
scan_fouriers(void)
{
int h, i, k;
k = tos;
i = 0;
while (fouriers[i]) {
h = tos;
while (fouriers[i]) {
scan(fouriers[i++]);
push_symbol(SYMBOL_A);
push(meta_a);
subst();
push_symbol(SYMBOL_B);
push(meta_b);
subst();
push_symbol(SYMBOL_C);
push(meta_c);
subst();
push_symbol(SYMBOL_N);
push(meta_n);
subst();
push_symbol(SYMBOL_X);
push(meta_x);
subst();
}
list(tos - h);
i++;
}
list(tos - k);
table_of_fourier = pop();
}
2007-05-08 16:57:30 +02:00
#if SELFTEST
2006-01-06 21:26:36 +01:00
static char *s[] = {
"#1",
2005-07-30 21:37:29 +02:00
"fourier(d(f(t,x),x),x)",
"i*x*fourier(f(t,x),x)",
2006-01-06 21:26:36 +01:00
"#2",
"eval(invfourier(i*x*fourier(f(t,x),x),x))",
"d(f(t,x),x)",
"#3",
2005-07-30 21:37:29 +02:00
"fourier(d(f(t,x),t),x)",
"d(fourier(f(t,x),x),t)",
2006-01-06 21:26:36 +01:00
"#4",
"eval(invfourier(d(fourier(f(t,x),x),t),x))",
"d(f(t,x),t)",
"#5",
2005-07-30 21:37:29 +02:00
"fourier(x*f(t,x),x)",
"i*d(fourier(f(t,x),x),x)",
2006-01-06 21:26:36 +01:00
"#6",
"eval(invfourier(i*d(fourier(f(t,x),x),x),x))",
"x*f(t,x)",
"#7",
"fourier(d(d(f(x),x),x))",
"-x^2*fourier(f(x),x)",
"#8",
"eval(invfourier(-x^2*fourier(f(x),x),x))",
"d(d(f(x),x),x)",
"#9",
"fourier(x^2*f(x),x)",
"-d(d(fourier(f(x),x),x),x)",
"#10",
"eval(invfourier(-d(d(fourier(f(x),x),x),x),x))",
"x^2*f(x)",
"#11",
2005-07-30 21:37:29 +02:00
"fourier(f(t,x)/x,x)",
"-i*integral(fourier(f(t,x),x),x)",
2006-01-06 21:26:36 +01:00
"#12",
"eval(invfourier(-i*integral(fourier(f(t,x),x),x),x))",
"f(t,x)/x",
"#13",
2005-07-30 21:37:29 +02:00
"fourier(integral(f(t,x),x),x)",
2006-01-06 21:26:36 +01:00
"-i*fourier(f(t,x),x)/x",
2005-07-30 21:37:29 +02:00
2006-01-06 21:26:36 +01:00
"#14",
"eval(invfourier(-i*fourier(f(t,x),x)/x,x))",
"integral(f(t,x),x)",
"#15",
"fourier(exp(-a^2*x^2),x)",
"pi^(1/2)*exp(-x^2/(4*a^2))/abs(a)",
"#16",
"eval(invfourier(pi^(1/2)*exp(-x^2/(4*a^2))/abs(a),x))",
"exp(-a^2*x^2)",
"#17",
"fourier(exp(-a*abs(x)),x)",
"2*a/(a^2+x^2)",
"#18",
"eval(invfourier(2*a/(a^2+x^2),x))",
"exp(-a*abs(x))",
"#19",
"fourier(heaviside(x)*exp(-a*x),x)",
"a/(a^2+x^2)-i*x/(a^2+x^2)",
"#20",
"eval(invfourier(a/(a^2+x^2)-i*x/(a^2+x^2),x))",
"1/2*exp(-a*abs(x))*sgn(x)+1/2*exp(-a*abs(x))",
2005-07-30 21:37:29 +02:00
2006-01-06 21:26:36 +01:00
"#21",
"fourier(sgn(x)*exp(-a*abs(x)),x)",
"-2*i*x/(a^2+x^2)",
"#22",
"eval(invfourier(-2*i*x/(a^2+x^2),x))",
"exp(-a*abs(x))*sgn(x)",
"#23",
"fourier(exp(i*y0*x-a*abs(x)),x)",
"2*a/(-2*x*y0+a^2+x^2+y0^2)",
"#24",
"fourier(heaviside(x)*exp(i*y0*x-a*x),x)",
"1/(a+i*x-i*y0)",
"#25",
"fourier(exp(i*x*y0)*carac(x,-L,L),x)",
"2*sin(L*x-L*y0)/(x-y0)",
"#26",
"fourier(sin(a*x),x)",
"i*pi*dirac(a+x)-i*pi*dirac(a-x)",
"#27",
"eval(invfourier(i*pi*dirac(a+x)-i*pi*dirac(a-x),x))",
"1/2*i*exp(-i*a*x)-1/2*i*exp(i*a*x)",
"#28",
"fourier(cos(a*x),x)",
"pi*dirac(a+x)+pi*dirac(a-x)",
"#29",
"eval(invfourier(pi*dirac(a+x)+pi*dirac(a-x),x))",
"1/2*exp(-i*a*x)+1/2*exp(i*a*x)",
"#30",
"fourier(sinh(a*x),x)",
"-pi*dirac(x-i*a)+pi*dirac(x+i*a)",
"#31",
"eval(invfourier(-pi*dirac(x-i*a)+pi*dirac(x+i*a),x))",
"-1/2*exp(-a*x)+1/2*exp(a*x)",
"#32",
"fourier(cosh(a*x),x)",
"pi*dirac(x-i*a)+pi*dirac(x+i*a)",
"#33",
"eval(invfourier(pi*dirac(x-i*a)+pi*dirac(x+i*a),x))",
"1/2*exp(-a*x)+1/2*exp(a*x)",
"#34",
"fourier(exp(i*b*x),x)",
"2*pi*dirac(b-x)",
"#35",
2005-07-30 21:37:29 +02:00
"fourier(dirac(x),x)",
"1",
2006-01-06 21:26:36 +01:00
"#36",
2005-07-30 21:37:29 +02:00
"fourier(1,x)",
"2*pi*dirac(x)",
2006-01-06 21:26:36 +01:00
"#37",
2005-07-30 21:37:29 +02:00
"fourier(d(dirac(A + x),x),x)",
"i*x*exp(i*A*x)",
2006-01-06 21:26:36 +01:00
"#38",
"fourier(carac(x,-1,1),x)",
"i*exp(-i*x)/x-i*exp(i*x)/x",
2005-07-30 21:37:29 +02:00
2006-01-06 21:26:36 +01:00
"#39",
"fourier(shah(a*x),x)",
"2*pi*shah(2*pi*x/a)/a",
2005-07-30 21:37:29 +02:00
2006-01-06 21:26:36 +01:00
"#40",
"fourier(1/x,x)",
"-i*pi*sgn(x)",
2005-07-30 21:37:29 +02:00
2006-01-06 21:26:36 +01:00
"#41",
"fourier(1/(a*x+b),x)",
"-i*pi*exp(i*b*x/a)*sgn(x)/a",
2005-07-30 21:37:29 +02:00
2006-01-06 21:26:36 +01:00
"#42",
"fourier(1/x^2,x)",
"-pi*x*sgn(x)",
2005-07-30 21:37:29 +02:00
2006-01-06 21:26:36 +01:00
"#43",
"fourier(1/(x^2+a^2),x)",
"pi*exp(-a*abs(x))/a",
2005-07-30 21:37:29 +02:00
2006-01-06 21:26:36 +01:00
"#44",
"fourier(1 / abs(x),x)",
"-2*euler-2*log(abs(x))",
2005-07-30 21:37:29 +02:00
2006-01-06 21:26:36 +01:00
"#45",
"fourier(sqrt(abs(x)))",
"-pi^(1/2)/(2^(1/2)*abs(x)^(3/2))",
2005-07-30 21:37:29 +02:00
2006-01-06 21:26:36 +01:00
"#46",
2005-07-30 21:37:29 +02:00
"fourier(exp(-a*abs(x)),x)",
"2*a/(a^2+x^2)",
2006-01-06 21:26:36 +01:00
"#47",
"eval(invfourier(2*a/(a^2+t^2),t))",
"exp(-a*abs(t))",
"#48",
2005-07-30 21:37:29 +02:00
"fourier(sgn(x)*exp(-a*abs(x)),x)",
"-2*i*x/(a^2+x^2)",
2006-01-06 21:26:36 +01:00
"#49",
"fourier(exp(-abs(x))/x,x)",
"-2*i*arctan(x)",
2005-07-30 21:37:29 +02:00
2006-01-06 21:26:36 +01:00
"#50",
"fourier(exp(-b abs(x) + i a x),x)",
"2*b/(-2*a*x+a^2+b^2+x^2)",
"#51",
"fourier(exp(-a x + i b x) heaviside(x),x)",
"1/(a-i*b+i*x)",
"#52",
"fourier(exp(-a*x^2),x)",
"pi^(1/2)*exp(-x^2/(4*a))/(abs(a)^(1/2))",
2005-10-04 02:02:16 +02:00
2006-01-06 21:26:36 +01:00
"#53",
"fourier(x*exp(-x^2))",
"-1/2*i*pi^(1/2)*x*exp(-1/4*x^2)",
2005-07-30 21:37:29 +02:00
2006-01-06 21:26:36 +01:00
"#54",
"rationalize(fourier(exp(i*a*x^2),x))",
"pi^(1/2)*exp(-i*x^2/(4*a))*(1+i*sgn(a))/(2^(1/2)*abs(a)^(1/2))",
"#55",
2005-07-30 21:37:29 +02:00
"fourier(besselj(x,0),x)",
"2/((1-x^2)^(1/2))",
2006-01-06 21:26:36 +01:00
"#56",
"eval(invfourier(2/((1-x^2)^(1/2)),x))",
"besselj(x,0)",
2005-07-30 21:37:29 +02:00
2006-01-06 21:26:36 +01:00
"#57",
"fourier(exp(i*b x) / (a^2 + x^2),x)",
"pi*exp(-a*abs(b-x))/a",
2005-07-30 21:37:29 +02:00
2006-01-06 21:26:36 +01:00
"#58",
"eval(fourier(dirac(t)-d(d(dirac(t),t),t),t))",
"1+t^2",
2005-07-30 21:37:29 +02:00
2006-01-06 21:26:36 +01:00
"#59",
"eval(fourier(fourier(d(dirac(t)*dirac(x),t)-d(d(dirac(t)*dirac(x),x),x),t),x))",
"i*t+x^2",
2005-07-30 21:37:29 +02:00
};
void
test_fourier(void)
{
test(__FILE__, s, sizeof s / sizeof (char *));
}
2007-05-08 16:57:30 +02:00
#endif