| Juhana Sadeharju on Fri, 31 Aug 2001 22:12:40 +0300 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
| plotminmax written/query for Jacobian functions |
Hello. I'm not in Pari-Gp development list but I follow the archives.
By the way, would you please set up digested versions of Pari/GP
mailing lists? They are useful if the lists are of only secondary
importance. While GP is only calculator I use frequently, I don't
have time to receive every mail separately.
-*-
Are there any plans to add Jacobi's elliptic functions and their
inverses to GP? I didn't found them. I could give a try if nobody
else has written them.
-*-
I wrote an alternative for plot() which could be named plotminmax().
The code is include at bottom, and it belongs to file "plotport.c".
I wrote plotminmax() because plot() most often only gives wrong curves.
The following example shows why I consider plotminmax() to be greatly
better -- but my implementation could be improved; see the code comments.
Magnitude of the frequency response of comb filter is given by h():
f(x,g,k) = x^(-k)/(1 - g*x^(-k));
h(t)=norm(f(exp(I*t),-0.7,10.0));
Which is then plotted with
plot(t=0.0,Pi,h(t))
and
plotminmax(t=0.0,Pi,h(t)).
plot() function:
10.900334 |''''''''''''''''''"''''''''''''''''''''''''"''''''''''''''''''|
| : : |
| : : |
| " : : " |
| : : : : |
| : : : : |
| : : __ : : |
| :: :: :: :: :: |
| :: : : :: : : :: |
| :: : : :: : : :: |
| : x : : :: : : x : |
| : : : : :: : : : : |
| : : : : : : : : : : |
| : : " : : : : " : : |
| : : : : : : : : : : |
| : : : " : : " : : : |
| x : : : : : : : : x |
| : : : : x x : : : : |
| : " _ : : _ " : |
| x " _ _ " x |
| _x "_ _" x_ __ __ _x "_ _" x_ |
0 "",,,,,,,,,"""",,,,,,,,,""",,,,,,,,,,""",,,,,,,,,"""",,,,,,,,,""
0 3.141593
plotminmax() function:
11.111111 |'''''*''''''''''''*'''''''''''**'''''''''''*''''''''''''*'''''|
| * ** ** ** * |
| * ** ** ** * |
| * ** ** ** * |
| * ** ** ** * |
| ** ** ** ** ** |
| ** ** ** ** ** |
| ** ** ** ** ** |
| *** ** ** ** *** |
| * * ** ** ** * * |
| * * ** ** ** * * |
| * * *** ** *** * * |
| * * * * ** * * * * |
| * * * * ** * * * * |
| * * * * **** * * * * |
| * * ** * * * * ** * * |
| ** * * * * * * * * ** |
| * ** * * * * * * ** * |
| * * * ** ** ** ** * * * |
| ** ** ** ** * * ** ** ** ** |
**** ******** ******** ******** ******** ****
0 |..............................................................|
0 3.141593
The plotminmax() shows both where the peaks are located and what are
their peak values, but plot() shows only the locations. plot() also
gives a misleading idea about the peak values, and that is not a good
thing at all.
So, finally you have someone who actually do use the crude plot()
function. ;-)
Best regards,
Juhana Sadeharju
PS. I'm that Juhana Kouhia in authors file; I implemented the initial
64-bit support. You could change the name, or at least remove the
invalid e-mail address.
==code for plotminmax()==
/*
In plotminmax() the function is heavily oversampled with respect to
screen resolution, and then minmax vertical lines are actually drawn
to screen. This plot style is more practical than the one used in
plot() which renders pretty much unusable due aliasing errors.
I render only with '*' characters now but plot()-style renderition
would work too. However, if we would plot to an image file, then we
could use an image-to-ascii converter program to downscale the image
to ascii text format with antialiasing.
*/
/* Oversampling ratio; it could be given as parameter to plotminmax()
* because 64 could be sometimes too much and sometimes too less.
*/
#define OSRATIO 64
/* fill_vertical() draws a vertical bar between min and max values
* within one column.
*/
static void
fill_vertical(screen scr, long i, int osmin, int osmax)
{
int j;
for(j = osmin; j <= osmax; j++) scr[i][j] = '*';
}
void
plotminmax(entree *ep, GEN a, GEN b, char *ch,GEN ysmlu,GEN ybigu, long prec)
{
long av = avma, av2,limite,j,i,sig;
GEN p1,p2,ysml,ybig,x,diff,dyj,dx,*y;
screen scr;
char buf[80];
int k,osmin,osmax;
y = (GEN *)malloc((OSRATIO*ISCR+2)*sizeof(GEN)); /* 2 includes one extra
sample */
sig=gcmp(b,a); if (!sig) return;
if (sig<0) { x=a; a=b; b=x; }
x = cgetr(prec); gaffect(a,x); push_val(ep, x);
for (i=1; i<=OSRATIO*ISCR+1; i++) y[i]=cgetr(3);
p1 = gdivgs(gsub(b,a), OSRATIO*ISCR); /* change done due one extra sample */
dx = cgetr(prec); gaffect(p1, dx);
/* NOTE: setting ysml and ybig this way leaves zero axis visible,
* and user is not able to use true min/max values for drawing range;
* it could be better to set ysml = +infty and ybig = -infty
* because then user is always able to set the other limit to
* zero;
* I have not done the change now. I wait for comments.
*/
ysml=gzero; ybig=gzero;
for (j=1; j<=JSCR; j++) scr[1][j]=scr[ISCR][j]=YY;
for (i=2; i<ISCR; i++)
{
scr[i][1] = XX_LOWER;
scr[i][JSCR]=XX_UPPER;
for (j=2; j<JSCR; j++) scr[i][j]=BLANK;
}
av2=avma; limite=stack_lim(av2,1);
/* NOTE: y[1] = f(a) and y[OSRATIO*ISCR+1] = f(b) */
for (i=1; i<=OSRATIO*ISCR+1; i++)
{
gaffect(READ_EXPR(ch),y[i]);
if (gcmp(y[i],ysml)<0) ysml=y[i];
if (gcmp(y[i],ybig)>0) ybig=y[i];
x = addrr(x,dx);
if (low_stack(limite, stack_lim(av2,1)))
{
long tetpil=avma;
if (DEBUGMEM>1) err(warnmem,"plot");
x = gerepile(av2,tetpil,rcopy(x));
}
ep->value = (void*)x;
}
if (ysmlu) ysml=ysmlu;
if (ybigu) ybig=ybigu;
avma=av2; diff=gsub(ybig,ysml);
if (gcmp0(diff)) { ybig=gaddsg(1,ybig); diff=gun; }
dyj = gdivsg(JSCR-1,diff);
av2=avma;
for (i=1; i<=ISCR; i++)
{
/* first character column spans values y[1] to y[OSRATIO+1];
* second character column spans values y[OSRATIO+1] to
* y[2*OSRATIO+1];
* last character column ISCR spans values y[(ISCR-1)*OSRATIO+1]
* to y[ISCR*OSRATIO+1]
*
* That is, there is one sample overlap between columns because
* the drawing should be continuous.
*/
osmin = 1000; /* unit is a screen character */
osmax = -1000;
for(k=(i-1)*OSRATIO+1; k<=i*OSRATIO+1; k++)
{
j = 1+gtolong(gmul(gsub(y[k],ysml),dyj));
if (osmin > j) osmin = j;
if (osmax < j) osmax = j;
}
fill_vertical(scr, i, osmin, osmax);
avma = av2;
}
p1=cgetr(3); gaffect(ybig,p1); pariputc('\n');
pariputsf("%s ", dsprintf9(rtodbl(p1), buf));
for (i=1; i<=ISCR; i++) pariputc(scr[i][JSCR]);
pariputc('\n');
for (j=(JSCR-1); j>=2; j--)
{
pariputs(" ");
for (i=1; i<=ISCR; i++) pariputc(scr[i][j]);
pariputc('\n');
}
p1=cgetr(3); gaffect(ysml,p1);
pariputsf("%s ", dsprintf9(rtodbl(p1), buf));
for (i=1; i<=ISCR; i++) pariputc(scr[i][1]);
pariputc('\n');
p1=cgetr(3); gaffect(a,p1);
p2=cgetr(3); gaffect(b,p2);
pariputsf("%10s%-9.7g%*.7g\n"," ",rtodbl(p1),ISCR-9,rtodbl(p2));
pop_val(ep); avma=av;
free(y);
}
==end==