xref: /plan9-contrib/sys/src/cmd/gs/src/zmath.c (revision 593dc095aefb2a85c828727bbfa9da139a49bdf4)
13ff48bf5SDavid du Colombier /* Copyright (C) 1989, 2000 Aladdin Enterprises.  All rights reserved.
27dd7cddfSDavid du Colombier 
3*593dc095SDavid du Colombier   This software is provided AS-IS with no warranty, either express or
4*593dc095SDavid du Colombier   implied.
57dd7cddfSDavid du Colombier 
6*593dc095SDavid du Colombier   This software is distributed under license and may not be copied,
7*593dc095SDavid du Colombier   modified or distributed except as expressly authorized under the terms
8*593dc095SDavid du Colombier   of the license contained in the file LICENSE in this distribution.
97dd7cddfSDavid du Colombier 
10*593dc095SDavid du Colombier   For more information about licensing, please refer to
11*593dc095SDavid du Colombier   http://www.ghostscript.com/licensing/. For information on
12*593dc095SDavid du Colombier   commercial licensing, go to http://www.artifex.com/licensing/ or
13*593dc095SDavid du Colombier   contact Artifex Software, Inc., 101 Lucas Valley Road #110,
14*593dc095SDavid du Colombier   San Rafael, CA  94903, U.S.A., +1(415)492-9861.
157dd7cddfSDavid du Colombier */
167dd7cddfSDavid du Colombier 
17*593dc095SDavid du Colombier /* $Id: zmath.c,v 1.5 2002/02/21 22:24:54 giles Exp $ */
187dd7cddfSDavid du Colombier /* Mathematical operators */
197dd7cddfSDavid du Colombier #include "math_.h"
207dd7cddfSDavid du Colombier #include "ghost.h"
217dd7cddfSDavid du Colombier #include "gxfarith.h"
227dd7cddfSDavid du Colombier #include "oper.h"
237dd7cddfSDavid du Colombier #include "store.h"
247dd7cddfSDavid du Colombier 
257dd7cddfSDavid du Colombier /*
267dd7cddfSDavid du Colombier  * Many of the procedures in this file are public only so they can be
277dd7cddfSDavid du Colombier  * called from the FunctionType 4 interpreter (zfunc4.c).
287dd7cddfSDavid du Colombier  */
297dd7cddfSDavid du Colombier 
307dd7cddfSDavid du Colombier /*
317dd7cddfSDavid du Colombier  * Define the current state of random number generator for operators.  We
327dd7cddfSDavid du Colombier  * have to implement this ourselves because the Unix rand doesn't provide
337dd7cddfSDavid du Colombier  * anything equivalent to rrand.  Note that the value always lies in the
347dd7cddfSDavid du Colombier  * range [0..0x7ffffffe], even if longs are longer than 32 bits.
357dd7cddfSDavid du Colombier  *
367dd7cddfSDavid du Colombier  * The state must be public so that context switching can save and
377dd7cddfSDavid du Colombier  * restore it.  (Even though the Red Book doesn't mention this,
387dd7cddfSDavid du Colombier  * we verified with Adobe that this is the case.)
397dd7cddfSDavid du Colombier  */
407dd7cddfSDavid du Colombier #define zrand_state (i_ctx_p->rand_state)
417dd7cddfSDavid du Colombier 
427dd7cddfSDavid du Colombier /* Initialize the random number generator. */
437dd7cddfSDavid du Colombier const long rand_state_initial = 1;
447dd7cddfSDavid du Colombier 
457dd7cddfSDavid du Colombier /****** NOTE: none of these operators currently ******/
467dd7cddfSDavid du Colombier /****** check for floating over- or underflow.	******/
477dd7cddfSDavid du Colombier 
487dd7cddfSDavid du Colombier /* <num> sqrt <real> */
497dd7cddfSDavid du Colombier int
zsqrt(i_ctx_t * i_ctx_p)507dd7cddfSDavid du Colombier zsqrt(i_ctx_t *i_ctx_p)
517dd7cddfSDavid du Colombier {
527dd7cddfSDavid du Colombier     os_ptr op = osp;
537dd7cddfSDavid du Colombier     double num;
547dd7cddfSDavid du Colombier     int code = real_param(op, &num);
557dd7cddfSDavid du Colombier 
567dd7cddfSDavid du Colombier     if (code < 0)
577dd7cddfSDavid du Colombier 	return code;
587dd7cddfSDavid du Colombier     if (num < 0.0)
597dd7cddfSDavid du Colombier 	return_error(e_rangecheck);
607dd7cddfSDavid du Colombier     make_real(op, sqrt(num));
617dd7cddfSDavid du Colombier     return 0;
627dd7cddfSDavid du Colombier }
637dd7cddfSDavid du Colombier 
647dd7cddfSDavid du Colombier /* <num> arccos <real> */
657dd7cddfSDavid du Colombier private int
zarccos(i_ctx_t * i_ctx_p)667dd7cddfSDavid du Colombier zarccos(i_ctx_t *i_ctx_p)
677dd7cddfSDavid du Colombier {
687dd7cddfSDavid du Colombier     os_ptr op = osp;
697dd7cddfSDavid du Colombier     double num, result;
707dd7cddfSDavid du Colombier     int code = real_param(op, &num);
717dd7cddfSDavid du Colombier 
727dd7cddfSDavid du Colombier     if (code < 0)
737dd7cddfSDavid du Colombier 	return code;
747dd7cddfSDavid du Colombier     result = acos(num) * radians_to_degrees;
757dd7cddfSDavid du Colombier     make_real(op, result);
767dd7cddfSDavid du Colombier     return 0;
777dd7cddfSDavid du Colombier }
787dd7cddfSDavid du Colombier 
797dd7cddfSDavid du Colombier /* <num> arcsin <real> */
807dd7cddfSDavid du Colombier private int
zarcsin(i_ctx_t * i_ctx_p)817dd7cddfSDavid du Colombier zarcsin(i_ctx_t *i_ctx_p)
827dd7cddfSDavid du Colombier {
837dd7cddfSDavid du Colombier     os_ptr op = osp;
847dd7cddfSDavid du Colombier     double num, result;
857dd7cddfSDavid du Colombier     int code = real_param(op, &num);
867dd7cddfSDavid du Colombier 
877dd7cddfSDavid du Colombier     if (code < 0)
887dd7cddfSDavid du Colombier 	return code;
897dd7cddfSDavid du Colombier     result = asin(num) * radians_to_degrees;
907dd7cddfSDavid du Colombier     make_real(op, result);
917dd7cddfSDavid du Colombier     return 0;
927dd7cddfSDavid du Colombier }
937dd7cddfSDavid du Colombier 
947dd7cddfSDavid du Colombier /* <num> <denom> atan <real> */
957dd7cddfSDavid du Colombier int
zatan(i_ctx_t * i_ctx_p)967dd7cddfSDavid du Colombier zatan(i_ctx_t *i_ctx_p)
977dd7cddfSDavid du Colombier {
987dd7cddfSDavid du Colombier     os_ptr op = osp;
997dd7cddfSDavid du Colombier     double args[2];
1007dd7cddfSDavid du Colombier     double result;
1017dd7cddfSDavid du Colombier     int code = num_params(op, 2, args);
1027dd7cddfSDavid du Colombier 
1037dd7cddfSDavid du Colombier     if (code < 0)
1047dd7cddfSDavid du Colombier 	return code;
1053ff48bf5SDavid du Colombier     code = gs_atan2_degrees(args[0], args[1], &result);
1063ff48bf5SDavid du Colombier     if (code < 0)
1073ff48bf5SDavid du Colombier 	return code;
1087dd7cddfSDavid du Colombier     make_real(op - 1, result);
1097dd7cddfSDavid du Colombier     pop(1);
1107dd7cddfSDavid du Colombier     return 0;
1117dd7cddfSDavid du Colombier }
1127dd7cddfSDavid du Colombier 
1137dd7cddfSDavid du Colombier /* <num> cos <real> */
1147dd7cddfSDavid du Colombier int
zcos(i_ctx_t * i_ctx_p)1157dd7cddfSDavid du Colombier zcos(i_ctx_t *i_ctx_p)
1167dd7cddfSDavid du Colombier {
1177dd7cddfSDavid du Colombier     os_ptr op = osp;
1187dd7cddfSDavid du Colombier     double angle;
1197dd7cddfSDavid du Colombier     int code = real_param(op, &angle);
1207dd7cddfSDavid du Colombier 
1217dd7cddfSDavid du Colombier     if (code < 0)
1227dd7cddfSDavid du Colombier 	return code;
1237dd7cddfSDavid du Colombier     make_real(op, gs_cos_degrees(angle));
1247dd7cddfSDavid du Colombier     return 0;
1257dd7cddfSDavid du Colombier }
1267dd7cddfSDavid du Colombier 
1277dd7cddfSDavid du Colombier /* <num> sin <real> */
1287dd7cddfSDavid du Colombier int
zsin(i_ctx_t * i_ctx_p)1297dd7cddfSDavid du Colombier zsin(i_ctx_t *i_ctx_p)
1307dd7cddfSDavid du Colombier {
1317dd7cddfSDavid du Colombier     os_ptr op = osp;
1327dd7cddfSDavid du Colombier     double angle;
1337dd7cddfSDavid du Colombier     int code = real_param(op, &angle);
1347dd7cddfSDavid du Colombier 
1357dd7cddfSDavid du Colombier     if (code < 0)
1367dd7cddfSDavid du Colombier 	return code;
1377dd7cddfSDavid du Colombier     make_real(op, gs_sin_degrees(angle));
1387dd7cddfSDavid du Colombier     return 0;
1397dd7cddfSDavid du Colombier }
1407dd7cddfSDavid du Colombier 
1417dd7cddfSDavid du Colombier /* <base> <exponent> exp <real> */
1427dd7cddfSDavid du Colombier int
zexp(i_ctx_t * i_ctx_p)1437dd7cddfSDavid du Colombier zexp(i_ctx_t *i_ctx_p)
1447dd7cddfSDavid du Colombier {
1457dd7cddfSDavid du Colombier     os_ptr op = osp;
1467dd7cddfSDavid du Colombier     double args[2];
1477dd7cddfSDavid du Colombier     double result;
1487dd7cddfSDavid du Colombier     double ipart;
1497dd7cddfSDavid du Colombier     int code = num_params(op, 2, args);
1507dd7cddfSDavid du Colombier 
1517dd7cddfSDavid du Colombier     if (code < 0)
1527dd7cddfSDavid du Colombier 	return code;
1537dd7cddfSDavid du Colombier     if (args[0] == 0.0 && args[1] == 0.0)
1547dd7cddfSDavid du Colombier 	return_error(e_undefinedresult);
1557dd7cddfSDavid du Colombier     if (args[0] < 0.0 && modf(args[1], &ipart) != 0.0)
1567dd7cddfSDavid du Colombier 	return_error(e_undefinedresult);
1577dd7cddfSDavid du Colombier     result = pow(args[0], args[1]);
1587dd7cddfSDavid du Colombier     make_real(op - 1, result);
1597dd7cddfSDavid du Colombier     pop(1);
1607dd7cddfSDavid du Colombier     return 0;
1617dd7cddfSDavid du Colombier }
1627dd7cddfSDavid du Colombier 
1637dd7cddfSDavid du Colombier /* <posnum> ln <real> */
1647dd7cddfSDavid du Colombier int
zln(i_ctx_t * i_ctx_p)1657dd7cddfSDavid du Colombier zln(i_ctx_t *i_ctx_p)
1667dd7cddfSDavid du Colombier {
1677dd7cddfSDavid du Colombier     os_ptr op = osp;
1687dd7cddfSDavid du Colombier     double num;
1697dd7cddfSDavid du Colombier     int code = real_param(op, &num);
1707dd7cddfSDavid du Colombier 
1717dd7cddfSDavid du Colombier     if (code < 0)
1727dd7cddfSDavid du Colombier 	return code;
1737dd7cddfSDavid du Colombier     if (num <= 0.0)
1747dd7cddfSDavid du Colombier 	return_error(e_rangecheck);
1757dd7cddfSDavid du Colombier     make_real(op, log(num));
1767dd7cddfSDavid du Colombier     return 0;
1777dd7cddfSDavid du Colombier }
1787dd7cddfSDavid du Colombier 
1797dd7cddfSDavid du Colombier /* <posnum> log <real> */
1807dd7cddfSDavid du Colombier int
zlog(i_ctx_t * i_ctx_p)1817dd7cddfSDavid du Colombier zlog(i_ctx_t *i_ctx_p)
1827dd7cddfSDavid du Colombier {
1837dd7cddfSDavid du Colombier     os_ptr op = osp;
1847dd7cddfSDavid du Colombier     double num;
1857dd7cddfSDavid du Colombier     int code = real_param(op, &num);
1867dd7cddfSDavid du Colombier 
1877dd7cddfSDavid du Colombier     if (code < 0)
1887dd7cddfSDavid du Colombier 	return code;
1897dd7cddfSDavid du Colombier     if (num <= 0.0)
1907dd7cddfSDavid du Colombier 	return_error(e_rangecheck);
1917dd7cddfSDavid du Colombier     make_real(op, log10(num));
1927dd7cddfSDavid du Colombier     return 0;
1937dd7cddfSDavid du Colombier }
1947dd7cddfSDavid du Colombier 
1957dd7cddfSDavid du Colombier /* - rand <int> */
1967dd7cddfSDavid du Colombier private int
zrand(i_ctx_t * i_ctx_p)1977dd7cddfSDavid du Colombier zrand(i_ctx_t *i_ctx_p)
1987dd7cddfSDavid du Colombier {
1997dd7cddfSDavid du Colombier     os_ptr op = osp;
2007dd7cddfSDavid du Colombier 
2017dd7cddfSDavid du Colombier 	/*
2027dd7cddfSDavid du Colombier 	 * We use an algorithm from CACM 31 no. 10, pp. 1192-1201,
2037dd7cddfSDavid du Colombier 	 * October 1988.  According to a posting by Ed Taft on
2047dd7cddfSDavid du Colombier 	 * comp.lang.postscript, Level 2 (Adobe) PostScript interpreters
2057dd7cddfSDavid du Colombier 	 * use this algorithm too:
2067dd7cddfSDavid du Colombier 	 *      x[n+1] = (16807 * x[n]) mod (2^31 - 1)
2077dd7cddfSDavid du Colombier 	 */
2087dd7cddfSDavid du Colombier #define A 16807
2097dd7cddfSDavid du Colombier #define M 0x7fffffff
2107dd7cddfSDavid du Colombier #define Q 127773		/* M / A */
2117dd7cddfSDavid du Colombier #define R 2836			/* M % A */
2127dd7cddfSDavid du Colombier     zrand_state = A * (zrand_state % Q) - R * (zrand_state / Q);
2137dd7cddfSDavid du Colombier     /* Note that zrand_state cannot be 0 here. */
2147dd7cddfSDavid du Colombier     if (zrand_state <= 0)
2157dd7cddfSDavid du Colombier 	zrand_state += M;
2167dd7cddfSDavid du Colombier #undef A
2177dd7cddfSDavid du Colombier #undef M
2187dd7cddfSDavid du Colombier #undef Q
2197dd7cddfSDavid du Colombier #undef R
2207dd7cddfSDavid du Colombier     push(1);
2217dd7cddfSDavid du Colombier     make_int(op, zrand_state);
2227dd7cddfSDavid du Colombier     return 0;
2237dd7cddfSDavid du Colombier }
2247dd7cddfSDavid du Colombier 
2257dd7cddfSDavid du Colombier /* <int> srand - */
2267dd7cddfSDavid du Colombier private int
zsrand(i_ctx_t * i_ctx_p)2277dd7cddfSDavid du Colombier zsrand(i_ctx_t *i_ctx_p)
2287dd7cddfSDavid du Colombier {
2297dd7cddfSDavid du Colombier     os_ptr op = osp;
2307dd7cddfSDavid du Colombier     long state;
2317dd7cddfSDavid du Colombier 
2327dd7cddfSDavid du Colombier     check_type(*op, t_integer);
2337dd7cddfSDavid du Colombier     state = op->value.intval;
2347dd7cddfSDavid du Colombier #if arch_sizeof_long > 4
2357dd7cddfSDavid du Colombier     /* Trim the state back to 32 bits. */
2367dd7cddfSDavid du Colombier     state = (int)state;
2377dd7cddfSDavid du Colombier #endif
2387dd7cddfSDavid du Colombier     /*
2397dd7cddfSDavid du Colombier      * The following somewhat bizarre adjustments are according to
2407dd7cddfSDavid du Colombier      * public information from Adobe describing their implementation.
2417dd7cddfSDavid du Colombier      */
2427dd7cddfSDavid du Colombier     if (state < 1)
2437dd7cddfSDavid du Colombier 	state = -(state % 0x7ffffffe) + 1;
2447dd7cddfSDavid du Colombier     else if (state > 0x7ffffffe)
2457dd7cddfSDavid du Colombier 	state = 0x7ffffffe;
2467dd7cddfSDavid du Colombier     zrand_state = state;
2477dd7cddfSDavid du Colombier     pop(1);
2487dd7cddfSDavid du Colombier     return 0;
2497dd7cddfSDavid du Colombier }
2507dd7cddfSDavid du Colombier 
2517dd7cddfSDavid du Colombier /* - rrand <int> */
2527dd7cddfSDavid du Colombier private int
zrrand(i_ctx_t * i_ctx_p)2537dd7cddfSDavid du Colombier zrrand(i_ctx_t *i_ctx_p)
2547dd7cddfSDavid du Colombier {
2557dd7cddfSDavid du Colombier     os_ptr op = osp;
2567dd7cddfSDavid du Colombier 
2577dd7cddfSDavid du Colombier     push(1);
2587dd7cddfSDavid du Colombier     make_int(op, zrand_state);
2597dd7cddfSDavid du Colombier     return 0;
2607dd7cddfSDavid du Colombier }
2617dd7cddfSDavid du Colombier 
2627dd7cddfSDavid du Colombier /* ------ Initialization procedure ------ */
2637dd7cddfSDavid du Colombier 
2647dd7cddfSDavid du Colombier const op_def zmath_op_defs[] =
2657dd7cddfSDavid du Colombier {
2667dd7cddfSDavid du Colombier     {"1arccos", zarccos},	/* extension */
2677dd7cddfSDavid du Colombier     {"1arcsin", zarcsin},	/* extension */
2687dd7cddfSDavid du Colombier     {"2atan", zatan},
2697dd7cddfSDavid du Colombier     {"1cos", zcos},
2707dd7cddfSDavid du Colombier     {"2exp", zexp},
2717dd7cddfSDavid du Colombier     {"1ln", zln},
2727dd7cddfSDavid du Colombier     {"1log", zlog},
2737dd7cddfSDavid du Colombier     {"0rand", zrand},
2747dd7cddfSDavid du Colombier     {"0rrand", zrrand},
2757dd7cddfSDavid du Colombier     {"1sin", zsin},
2767dd7cddfSDavid du Colombier     {"1sqrt", zsqrt},
2777dd7cddfSDavid du Colombier     {"1srand", zsrand},
2787dd7cddfSDavid du Colombier     op_def_end(0)
2797dd7cddfSDavid du Colombier };
280