1 /* Copyright (C) 2001-2019 Artifex Software, Inc.
2    All Rights Reserved.
3 
4    This software is provided AS-IS with no warranty, either express or
5    implied.
6 
7    This software is distributed under license and may not be copied,
8    modified or distributed except as expressly authorized under the terms
9    of the license contained in the file LICENSE in this distribution.
10 
11    Refer to licensing information at http://www.artifex.com or contact
12    Artifex Software, Inc.,  1305 Grant Avenue - Suite 200, Novato,
13    CA 94945, U.S.A., +1(415)492-9861, for further information.
14 */
15 
16 
17 /* Double-precision floating point arithmetic operators */
18 #include "math_.h"
19 #include "memory_.h"
20 #include "string_.h"
21 #include "ctype_.h"
22 #include "ghost.h"
23 #include "gxfarith.h"
24 #include "oper.h"
25 #include "store.h"
26 
27 /*
28  * Thanks to Jean-Pierre Demailly of the Institut Fourier of the
29  * Universit\'e de Grenoble I <demailly@fourier.grenet.fr> for proposing
30  * this package and for arranging the funding for its creation.
31  *
32  * These operators work with doubles represented as 8-byte strings.  When
33  * applicable, they write their result into a string supplied as an argument.
34  * They also accept ints and reals as arguments.
35  */
36 
37 /* Forward references */
38 static int double_params_result(os_ptr, int, double *);
39 static int double_params(os_ptr, int, double *);
40 static int double_result(i_ctx_t *, int, double);
41 static int double_unary(i_ctx_t *, double (*)(double));
42 
43 #define dbegin_unary()\
44         os_ptr op = osp;\
45         double num;\
46         int code = double_params_result(op, 1, &num);\
47 \
48         if ( code < 0 )\
49           return code
50 
51 #define dbegin_binary()\
52         os_ptr op = osp;\
53         double num[2];\
54         int code = double_params_result(op, 2, num);\
55 \
56         if ( code < 0 )\
57           return code
58 
59 /* ------ Arithmetic ------ */
60 
61 /* <dnum1> <dnum2> <dresult> .dadd <dresult> */
62 static int
zdadd(i_ctx_t * i_ctx_p)63 zdadd(i_ctx_t *i_ctx_p)
64 {
65     dbegin_binary();
66     return double_result(i_ctx_p, 2, num[0] + num[1]);
67 }
68 
69 /* <dnum1> <dnum2> <dresult> .ddiv <dresult> */
70 static int
zddiv(i_ctx_t * i_ctx_p)71 zddiv(i_ctx_t *i_ctx_p)
72 {
73     dbegin_binary();
74     if (num[1] == 0.0)
75         return_error(gs_error_undefinedresult);
76     return double_result(i_ctx_p, 2, num[0] / num[1]);
77 }
78 
79 /* <dnum1> <dnum2> <dresult> .dmul <dresult> */
80 static int
zdmul(i_ctx_t * i_ctx_p)81 zdmul(i_ctx_t *i_ctx_p)
82 {
83     dbegin_binary();
84     return double_result(i_ctx_p, 2, num[0] * num[1]);
85 }
86 
87 /* <dnum1> <dnum2> <dresult> .dsub <dresult> */
88 static int
zdsub(i_ctx_t * i_ctx_p)89 zdsub(i_ctx_t *i_ctx_p)
90 {
91     dbegin_binary();
92     return double_result(i_ctx_p, 2, num[0] - num[1]);
93 }
94 
95 /* ------ Simple functions ------ */
96 
97 /* <dnum> <dresult> .dabs <dresult> */
98 static int
zdabs(i_ctx_t * i_ctx_p)99 zdabs(i_ctx_t *i_ctx_p)
100 {
101     return double_unary(i_ctx_p, fabs);
102 }
103 
104 /* <dnum> <dresult> .dceiling <dresult> */
105 static int
zdceiling(i_ctx_t * i_ctx_p)106 zdceiling(i_ctx_t *i_ctx_p)
107 {
108     return double_unary(i_ctx_p, ceil);
109 }
110 
111 /* <dnum> <dresult> .dfloor <dresult> */
112 static int
zdfloor(i_ctx_t * i_ctx_p)113 zdfloor(i_ctx_t *i_ctx_p)
114 {
115     return double_unary(i_ctx_p, floor);
116 }
117 
118 /* <dnum> <dresult> .dneg <dresult> */
119 static int
zdneg(i_ctx_t * i_ctx_p)120 zdneg(i_ctx_t *i_ctx_p)
121 {
122     dbegin_unary();
123     return double_result(i_ctx_p, 1, -num);
124 }
125 
126 /* <dnum> <dresult> .dround <dresult> */
127 static int
zdround(i_ctx_t * i_ctx_p)128 zdround(i_ctx_t *i_ctx_p)
129 {
130     dbegin_unary();
131     return double_result(i_ctx_p, 1, floor(num + 0.5));
132 }
133 
134 /* <dnum> <dresult> .dsqrt <dresult> */
135 static int
zdsqrt(i_ctx_t * i_ctx_p)136 zdsqrt(i_ctx_t *i_ctx_p)
137 {
138     dbegin_unary();
139     if (num < 0.0)
140         return_error(gs_error_rangecheck);
141     return double_result(i_ctx_p, 1, sqrt(num));
142 }
143 
144 /* <dnum> <dresult> .dtruncate <dresult> */
145 static int
zdtruncate(i_ctx_t * i_ctx_p)146 zdtruncate(i_ctx_t *i_ctx_p)
147 {
148     dbegin_unary();
149     return double_result(i_ctx_p, 1, (num < 0 ? ceil(num) : floor(num)));
150 }
151 
152 /* ------ Transcendental functions ------ */
153 
154 static int
darc(i_ctx_t * i_ctx_p,double (* afunc)(double))155 darc(i_ctx_t *i_ctx_p, double (*afunc)(double))
156 {
157     dbegin_unary();
158     return double_result(i_ctx_p, 1, (*afunc)(num) * radians_to_degrees);
159 }
160 /* <dnum> <dresult> .darccos <dresult> */
161 static int
zdarccos(i_ctx_t * i_ctx_p)162 zdarccos(i_ctx_t *i_ctx_p)
163 {
164     return darc(i_ctx_p, acos);
165 }
166 /* <dnum> <dresult> .darcsin <dresult> */
167 static int
zdarcsin(i_ctx_t * i_ctx_p)168 zdarcsin(i_ctx_t *i_ctx_p)
169 {
170     return darc(i_ctx_p, asin);
171 }
172 
173 /* <dnum> <ddenom> <dresult> .datan <dresult> */
174 static int
zdatan(i_ctx_t * i_ctx_p)175 zdatan(i_ctx_t *i_ctx_p)
176 {
177     double result;
178 
179     dbegin_binary();
180     if (num[0] == 0) {		/* on X-axis, special case */
181         if (num[1] == 0)
182             return_error(gs_error_undefinedresult);
183         result = (num[1] < 0 ? 180 : 0);
184     } else {
185         result = atan2(num[0], num[1]) * radians_to_degrees;
186         if (result < 0)
187             result += 360;
188     }
189     return double_result(i_ctx_p, 2, result);
190 }
191 
192 /* <dnum> <dresult> .dcos <dresult> */
193 static int
zdcos(i_ctx_t * i_ctx_p)194 zdcos(i_ctx_t *i_ctx_p)
195 {
196     return double_unary(i_ctx_p, gs_cos_degrees);
197 }
198 
199 /* <dbase> <dexponent> <dresult> .dexp <dresult> */
200 static int
zdexp(i_ctx_t * i_ctx_p)201 zdexp(i_ctx_t *i_ctx_p)
202 {
203     double ipart;
204 
205     dbegin_binary();
206     if (num[0] == 0.0 && num[1] == 0.0)
207         return_error(gs_error_undefinedresult);
208     if (num[0] < 0.0 && modf(num[1], &ipart) != 0.0)
209         return_error(gs_error_undefinedresult);
210     return double_result(i_ctx_p, 2, pow(num[0], num[1]));
211 }
212 
213 static int
dlog(i_ctx_t * i_ctx_p,double (* lfunc)(double))214 dlog(i_ctx_t *i_ctx_p, double (*lfunc)(double))
215 {
216     dbegin_unary();
217     if (num <= 0.0)
218         return_error(gs_error_rangecheck);
219     return double_result(i_ctx_p, 1, (*lfunc)(num));
220 }
221 /* <dposnum> <dresult> .dln <dresult> */
222 static int
zdln(i_ctx_t * i_ctx_p)223 zdln(i_ctx_t *i_ctx_p)
224 {
225     return dlog(i_ctx_p, log);
226 }
227 /* <dposnum> <dresult> .dlog <dresult> */
228 static int
zdlog(i_ctx_t * i_ctx_p)229 zdlog(i_ctx_t *i_ctx_p)
230 {
231     return dlog(i_ctx_p, log10);
232 }
233 
234 /* <dnum> <dresult> .dsin <dresult> */
235 static int
zdsin(i_ctx_t * i_ctx_p)236 zdsin(i_ctx_t *i_ctx_p)
237 {
238     return double_unary(i_ctx_p, gs_sin_degrees);
239 }
240 
241 /* ------ Comparison ------ */
242 
243 static int
dcompare(i_ctx_t * i_ctx_p,int mask)244 dcompare(i_ctx_t *i_ctx_p, int mask)
245 {
246     os_ptr op = osp;
247     double num[2];
248     int code = double_params(op, 2, num);
249 
250     if (code < 0)
251         return code;
252     make_bool(op - 1,
253               (mask & (num[0] < num[1] ? 1 : num[0] > num[1] ? 4 : 2))
254               != 0);
255     pop(1);
256     return 0;
257 }
258 /* <dnum1> <dnum2> .deq <bool> */
259 static int
zdeq(i_ctx_t * i_ctx_p)260 zdeq(i_ctx_t *i_ctx_p)
261 {
262     return dcompare(i_ctx_p, 2);
263 }
264 /* <dnum1> <dnum2> .dge <bool> */
265 static int
zdge(i_ctx_t * i_ctx_p)266 zdge(i_ctx_t *i_ctx_p)
267 {
268     return dcompare(i_ctx_p, 6);
269 }
270 /* <dnum1> <dnum2> .dgt <bool> */
271 static int
zdgt(i_ctx_t * i_ctx_p)272 zdgt(i_ctx_t *i_ctx_p)
273 {
274     return dcompare(i_ctx_p, 4);
275 }
276 /* <dnum1> <dnum2> .dle <bool> */
277 static int
zdle(i_ctx_t * i_ctx_p)278 zdle(i_ctx_t *i_ctx_p)
279 {
280     return dcompare(i_ctx_p, 3);
281 }
282 /* <dnum1> <dnum2> .dlt <bool> */
283 static int
zdlt(i_ctx_t * i_ctx_p)284 zdlt(i_ctx_t *i_ctx_p)
285 {
286     return dcompare(i_ctx_p, 1);
287 }
288 /* <dnum1> <dnum2> .dne <bool> */
289 static int
zdne(i_ctx_t * i_ctx_p)290 zdne(i_ctx_t *i_ctx_p)
291 {
292     return dcompare(i_ctx_p, 5);
293 }
294 
295 /* ------ Conversion ------ */
296 
297 /* Take the easy way out.... */
298 #define MAX_CHARS 50
299 
300 /* <dnum> <dresult> .cvd <dresult> */
301 static int
zcvd(i_ctx_t * i_ctx_p)302 zcvd(i_ctx_t *i_ctx_p)
303 {
304     dbegin_unary();
305     return double_result(i_ctx_p, 1, num);
306 }
307 
308 /* <string> <dresult> .cvsd <dresult> */
309 static int
zcvsd(i_ctx_t * i_ctx_p)310 zcvsd(i_ctx_t *i_ctx_p)
311 {
312     os_ptr op = osp;
313     int code = double_params_result(op, 0, NULL);
314     double num;
315     char dot, buf[MAX_CHARS + 2];
316     char *str = buf;
317     uint len;
318     char end;
319 
320     if (code < 0)
321         return code;
322     check_read_type(op[-1], t_string);
323     len = r_size(op - 1);
324     if (len > MAX_CHARS)
325         return_error(gs_error_limitcheck);
326     gs_sprintf(buf, "%f", 1.5);
327     dot = buf[1]; /* locale-dependent */
328     memcpy(str, op[-1].value.bytes, len);
329     /*
330      * We check syntax in the following way: we remove whitespace,
331      * verify that the string contains only [0123456789+-.dDeE],
332      * then append a $ and then check that the next character after
333      * the scanned number is a $.
334      */
335     while (len > 0 && isspace(*str))
336         ++str, --len;
337     while (len > 0 && isspace(str[len - 1]))
338         --len;
339     str[len] = 0;
340     if (strspn(str, "0123456789+-.dDeE") != len)
341         return_error(gs_error_syntaxerror);
342     strcat(str, "$");
343     if (dot != '.') {
344         char *pdot = strchr(str, '.');
345         if (pdot)
346             *pdot = dot;
347     }
348     if (sscanf(str, "%lf%c", &num, &end) != 2 || end != '$')
349         return_error(gs_error_syntaxerror);
350     return double_result(i_ctx_p, 1, num);
351 }
352 
353 /* <dnum> .dcvi <int> */
354 static int
zdcvi(i_ctx_t * i_ctx_p)355 zdcvi(i_ctx_t *i_ctx_p)
356 {
357     os_ptr op = osp;
358 #define alt_min_long (-1L << (ARCH_SIZEOF_LONG * 8 - 1))
359 #define alt_max_long (~(alt_min_long))
360     static const double min_int_real = (alt_min_long * 1.0 - 1);
361     static const double max_int_real = (alt_max_long * 1.0 + 1);
362     double num;
363     int code = double_params(op, 1, &num);
364 
365     if (code < 0)
366         return code;
367 
368     if (num < min_int_real || num > max_int_real)
369         return_error(gs_error_rangecheck);
370     make_int(op, (long)num);	/* truncates toward 0 */
371     return 0;
372 }
373 
374 /* <dnum> .dcvr <real> */
375 static int
zdcvr(i_ctx_t * i_ctx_p)376 zdcvr(i_ctx_t *i_ctx_p)
377 {
378     os_ptr op = osp;
379 #define b30 (0x40000000L * 1.0)
380 #define max_mag (0xffffff * b30 * b30 * b30 * 0x4000)
381     static const float min_real = -max_mag;
382     static const float max_real = max_mag;
383 #undef b30
384 #undef max_mag
385     double num;
386     int code = double_params(op, 1, &num);
387 
388     if (code < 0)
389         return code;
390     if (num < min_real || num > max_real)
391         return_error(gs_error_rangecheck);
392     make_real(op, (float)num);
393     return 0;
394 }
395 
396 /* <dnum> <string> .dcvs <substring> */
397 static int
zdcvs(i_ctx_t * i_ctx_p)398 zdcvs(i_ctx_t *i_ctx_p)
399 {
400     os_ptr op = osp;
401     double num;
402     int code = double_params(op - 1, 1, &num);
403     char dot, str[MAX_CHARS + 1];
404     int len;
405 
406     if (code < 0)
407         return code;
408     check_write_type(*op, t_string);
409     gs_sprintf(str, "%f", 1.5);
410     dot = str[1]; /* locale-dependent */
411     /*
412      * To get fully accurate output results for IEEE double-
413      * precision floats (53 bits of mantissa), the ANSI
414      * %g default of 6 digits is not enough; 16 are needed.
415      * Unfortunately, using %.16g produces unfortunate artifacts such as
416      * 1.2 printing as 1.200000000000005.  Therefore, we print using %g,
417      * and if the result isn't accurate enough, print again
418      * using %.16g.
419      */
420     {
421         double scanned;
422 
423         gs_sprintf(str, "%g", num);
424         sscanf(str, "%lf", &scanned);
425         if (scanned != num)
426             gs_sprintf(str, "%.16g", num);
427     }
428     len = strlen(str);
429     if (len > r_size(op))
430         return_error(gs_error_rangecheck);
431     /* Juggling locales isn't thread-safe. Posix me harder. */
432     if (dot != '.') {
433         char *pdot = strchr(str, dot);
434         if (pdot)
435             *pdot = '.';
436     }
437     memcpy(op->value.bytes, str, len);
438     op[-1] = *op;
439     r_set_size(op - 1, len);
440     pop(1);
441     return 0;
442 }
443 
444 /* ------ Initialization table ------ */
445 
446 /* We need to split the table because of the 16-element limit. */
447 const op_def zdouble1_op_defs[] = {
448                 /* Arithmetic */
449     {"3.dadd", zdadd},
450     {"3.ddiv", zddiv},
451     {"3.dmul", zdmul},
452     {"3.dsub", zdsub},
453                 /* Comparison */
454     {"2.deq", zdeq},
455     {"2.dge", zdge},
456     {"2.dgt", zdgt},
457     {"2.dle", zdle},
458     {"2.dlt", zdlt},
459     {"2.dne", zdne},
460                 /* Conversion */
461     {"2.cvd", zcvd},
462     {"2.cvsd", zcvsd},
463     {"1.dcvi", zdcvi},
464     {"1.dcvr", zdcvr},
465     {"2.dcvs", zdcvs},
466     op_def_end(0)
467 };
468 const op_def zdouble2_op_defs[] = {
469                 /* Simple functions */
470     {"2.dabs", zdabs},
471     {"2.dceiling", zdceiling},
472     {"2.dfloor", zdfloor},
473     {"2.dneg", zdneg},
474     {"2.dround", zdround},
475     {"2.dsqrt", zdsqrt},
476     {"2.dtruncate", zdtruncate},
477                 /* Transcendental functions */
478     {"2.darccos", zdarccos},
479     {"2.darcsin", zdarcsin},
480     {"3.datan", zdatan},
481     {"2.dcos", zdcos},
482     {"3.dexp", zdexp},
483     {"2.dln", zdln},
484     {"2.dlog", zdlog},
485     {"2.dsin", zdsin},
486     op_def_end(0)
487 };
488 
489 /* ------ Internal procedures ------ */
490 
491 /* Get some double arguments. */
492 static int
double_params(os_ptr op,int count,double * pval)493 double_params(os_ptr op, int count, double *pval)
494 {
495     pval += count;
496     while (--count >= 0) {
497         switch (r_type(op)) {
498             case t_real:
499                 *--pval = op->value.realval;
500                 break;
501             case t_integer:
502                 *--pval = op->value.intval;
503                 break;
504             case t_string:
505                 if (!r_has_attr(op, a_read) ||
506                     r_size(op) != sizeof(double)
507                 )
508                            return_error(gs_error_typecheck);
509                 --pval;
510                 memcpy(pval, op->value.bytes, sizeof(double));
511                 break;
512             case t__invalid:
513                 return_error(gs_error_stackunderflow);
514             default:
515                 return_error(gs_error_typecheck);
516         }
517         op--;
518     }
519     return 0;
520 }
521 
522 /* Get some double arguments, and check for a double result. */
523 static int
double_params_result(os_ptr op,int count,double * pval)524 double_params_result(os_ptr op, int count, double *pval)
525 {
526     check_write_type(*op, t_string);
527     if (r_size(op) != sizeof(double))
528         return_error(gs_error_typecheck);
529     return double_params(op - 1, count, pval);
530 }
531 
532 /* Return a double result. */
533 static int
double_result(i_ctx_t * i_ctx_p,int count,double result)534 double_result(i_ctx_t *i_ctx_p, int count, double result)
535 {
536     os_ptr op = osp;
537     os_ptr op1 = op - count;
538 
539     ref_assign_inline(op1, op);
540     memcpy(op1->value.bytes, &result, sizeof(double));
541     pop(count);
542     return 0;
543 }
544 
545 /* Apply a unary function to a double operand. */
546 static int
double_unary(i_ctx_t * i_ctx_p,double (* func)(double))547 double_unary(i_ctx_t *i_ctx_p, double (*func)(double))
548 {
549     dbegin_unary();
550     return double_result(i_ctx_p, 1, (*func)(num));
551 }
552