Line data Source code
1 : /* Copyright (C) 2016 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : /********************************************************************/
19 : /** **/
20 : /** LOW-RES PLOT **/
21 : /** **/
22 : /********************************************************************/
23 : #define ISCR 64
24 : #define JSCR 22
25 :
26 : INLINE long
27 1920 : DTOL(double t) { return (long)(t + 0.5); }
28 :
29 : static char
30 1920 : PICT(long j) {
31 1920 : switch(j%3) {
32 862 : case 0: return '_';
33 368 : case 1: return 'x';
34 690 : default: return '"';
35 : }
36 : }
37 : static char
38 30 : PICTZERO(long j) {
39 30 : switch(j%3) {
40 8 : case 0: return ',';
41 11 : case 1: return '-';
42 11 : default: return '`';
43 : }
44 : }
45 :
46 : static char *
47 60 : dsprintf9(double d, char *buf)
48 : {
49 60 : int i = 10;
50 :
51 115 : while (--i >= 0) {
52 115 : sprintf(buf, "%9.*g", i, d);
53 115 : if (strlen(buf) <= 9) break;
54 : }
55 60 : return buf;
56 : }
57 :
58 : typedef unsigned char screen[ISCR+1][JSCR+1];
59 :
60 : static void
61 1890 : fill_gap(screen scr, long i, int jnew, int jpre)
62 : {
63 : int mid, i_up, i_lo, up, lo;
64 :
65 1890 : if (jpre < jnew - 2) {
66 60 : up = jnew - 1; i_up = i;
67 60 : lo = jpre + 1; i_lo = i - 1;
68 1830 : } else if (jnew < jpre - 2) {
69 88 : up = jpre - 1; i_up = i - 1;
70 88 : lo = jnew + 1; i_lo = i;
71 1742 : } else return; /* if gap < 2, leave it as it is. */
72 :
73 148 : mid = (jpre+jnew)/2;
74 148 : if (mid>JSCR) mid=JSCR; else if (mid<0) mid=0;
75 148 : if (lo<0) lo=0;
76 438 : if (lo<=JSCR) while (lo <= mid) scr[i_lo][lo++] = ':';
77 148 : if (up>JSCR) up=JSCR;
78 424 : if (up>=0) while (up > mid) scr[i_up][up--] = ':';
79 : }
80 :
81 : static double
82 60 : todbl(GEN x) { return rtodbl(gtofp(x, LOWDEFAULTPREC)); }
83 :
84 : void
85 30 : pariplot(void* E, GEN (*fun)(void *E, GEN x), GEN a, GEN b, GEN ysmlu,GEN ybigu, long prec)
86 : {
87 30 : const char BLANK = ' ', YY = '|', XX_UPPER = '\'', XX_LOWER = '.';
88 : long jz, j, i, sig;
89 30 : pari_sp av = avma;
90 30 : int jnew, jpre = 0; /* for lint */
91 : GEN x, dx;
92 : double diff, dyj, ysml, ybig, y[ISCR+1];
93 : screen scr;
94 : char buf[80], z;
95 :
96 30 : sig=gcmp(b,a); if (!sig) return;
97 30 : if (sig<0) { x=a; a=b; b=x; }
98 30 : x = gtofp(a, prec);
99 30 : dx = divru(gtofp(gsub(b,a),prec), ISCR-1);
100 690 : for (j=1; j<=JSCR; j++) scr[1][j]=scr[ISCR][j]=YY;
101 1890 : for (i=2; i<ISCR; i++)
102 : {
103 1860 : scr[i][1] = XX_LOWER;
104 1860 : scr[i][JSCR]= XX_UPPER;
105 39060 : for (j=2; j<JSCR; j++) scr[i][j] = BLANK;
106 : }
107 30 : ysml = ybig = 0.; /* -Wall */
108 1950 : for (i=1; i<=ISCR; i++)
109 : {
110 1920 : pari_sp av2 = avma;
111 1920 : y[i] = gtodouble( fun(E, x) );
112 1920 : set_avma(av2);
113 1920 : if (i == 1)
114 30 : ysml = ybig = y[1];
115 : else
116 : {
117 1890 : if (y[i] < ysml) ysml = y[i];
118 1890 : if (y[i] > ybig) ybig = y[i];
119 : }
120 1920 : x = addrr(x,dx);
121 : }
122 30 : set_avma(av);
123 30 : if (ysmlu) ysml = gtodouble(ysmlu);
124 30 : if (ybigu) ybig = gtodouble(ybigu);
125 30 : diff = ybig - ysml;
126 30 : if (!diff) { ybig += 1; diff= 1.; }
127 30 : dyj = ((JSCR-1)*3+2) / diff;
128 : /* work around bug in gcc-4.8 (32bit): plot(x=-5,5,sin(x)))) */
129 30 : jz = 3 - (long)(ysml*dyj + 0.5); /* 3 - DTOL(ysml*dyj) */
130 30 : z = PICTZERO(jz); jz /= 3;
131 1950 : for (i=1; i<=ISCR; i++)
132 : {
133 1920 : if (0<=jz && jz<=JSCR) scr[i][jz]=z;
134 1920 : j = 3 + DTOL((y[i]-ysml)*dyj);
135 1920 : jnew = j/3;
136 1920 : if (i > 1) fill_gap(scr, i, jnew, jpre);
137 1920 : if (0<=jnew && jnew<=JSCR) scr[i][jnew] = PICT(j);
138 1920 : jpre = jnew;
139 : }
140 30 : pari_putc('\n');
141 30 : pari_printf("%s ", dsprintf9(ybig, buf));
142 1950 : for (i=1; i<=ISCR; i++) pari_putc(scr[i][JSCR]);
143 30 : pari_putc('\n');
144 630 : for (j=(JSCR-1); j>=2; j--)
145 : {
146 600 : pari_puts(" ");
147 39000 : for (i=1; i<=ISCR; i++) pari_putc(scr[i][j]);
148 600 : pari_putc('\n');
149 : }
150 30 : pari_printf("%s ", dsprintf9(ysml, buf));
151 1950 : for (i=1; i<=ISCR; i++) pari_putc(scr[i][1]);
152 30 : pari_putc('\n');
153 : {
154 : char line[10 + 32 + 32 + ISCR - 9];
155 30 : sprintf(line, "%10s%-9.7g%*.7g\n"," ",todbl(a),ISCR-9,todbl(b));
156 30 : pari_printf(line);
157 : }
158 : }
159 :
160 : void
161 30 : pariplot0(GEN a, GEN b, GEN code, GEN ysmlu,GEN ybigu, long prec)
162 : {
163 30 : push_lex(gen_0, code);
164 30 : pariplot((void*)code, &gp_eval, a, b, ysmlu, ybigu, prec);
165 30 : pop_lex(1);
166 30 : }
|