14a238c70SJohn Marino /* mpfr_exp -- exponential of a floating-point number
24a238c70SJohn Marino
3*ab6d115fSJohn Marino Copyright 1999, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4*ab6d115fSJohn Marino Contributed by the AriC and Caramel projects, INRIA.
54a238c70SJohn Marino
64a238c70SJohn Marino This file is part of the GNU MPFR Library.
74a238c70SJohn Marino
84a238c70SJohn Marino The GNU MPFR Library is free software; you can redistribute it and/or modify
94a238c70SJohn Marino it under the terms of the GNU Lesser General Public License as published by
104a238c70SJohn Marino the Free Software Foundation; either version 3 of the License, or (at your
114a238c70SJohn Marino option) any later version.
124a238c70SJohn Marino
134a238c70SJohn Marino The GNU MPFR Library is distributed in the hope that it will be useful, but
144a238c70SJohn Marino WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
154a238c70SJohn Marino or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
164a238c70SJohn Marino License for more details.
174a238c70SJohn Marino
184a238c70SJohn Marino You should have received a copy of the GNU Lesser General Public License
194a238c70SJohn Marino along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
204a238c70SJohn Marino http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
214a238c70SJohn Marino 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
224a238c70SJohn Marino
234a238c70SJohn Marino #define MPFR_NEED_LONGLONG_H /* for MPFR_MPZ_SIZEINBASE2 */
244a238c70SJohn Marino #include "mpfr-impl.h"
254a238c70SJohn Marino
264a238c70SJohn Marino /* y <- exp(p/2^r) within 1 ulp, using 2^m terms from the series
274a238c70SJohn Marino Assume |p/2^r| < 1.
284a238c70SJohn Marino We use the following binary splitting formula:
294a238c70SJohn Marino P(a,b) = p if a+1=b, P(a,c)*P(c,b) otherwise
304a238c70SJohn Marino Q(a,b) = a*2^r if a+1=b [except Q(0,1)=1], Q(a,c)*Q(c,b) otherwise
314a238c70SJohn Marino T(a,b) = P(a,b) if a+1=b, Q(c,b)*T(a,c)+P(a,c)*T(c,b) otherwise
324a238c70SJohn Marino Then exp(p/2^r) ~ T(0,i)/Q(0,i) for i so that (p/2^r)^i/i! is small enough.
334a238c70SJohn Marino
344a238c70SJohn Marino Since P(a,b) = p^(b-a), and we consider only values of b-a of the form 2^j,
354a238c70SJohn Marino we don't need to compute P(), we only precompute p^(2^j) in the ptoj[] array
364a238c70SJohn Marino below.
374a238c70SJohn Marino
384a238c70SJohn Marino Since Q(a,b) is divisible by 2^(r*(b-a-1)), we don't compute the power of
394a238c70SJohn Marino two part.
404a238c70SJohn Marino */
414a238c70SJohn Marino static void
mpfr_exp_rational(mpfr_ptr y,mpz_ptr p,long r,int m,mpz_t * Q,mpfr_prec_t * mult)424a238c70SJohn Marino mpfr_exp_rational (mpfr_ptr y, mpz_ptr p, long r, int m,
434a238c70SJohn Marino mpz_t *Q, mpfr_prec_t *mult)
444a238c70SJohn Marino {
454a238c70SJohn Marino unsigned long n, i, j;
464a238c70SJohn Marino mpz_t *S, *ptoj;
474a238c70SJohn Marino mpfr_prec_t *log2_nb_terms;
484a238c70SJohn Marino mpfr_exp_t diff, expo;
494a238c70SJohn Marino mpfr_prec_t precy = MPFR_PREC(y), prec_i_have, prec_ptoj;
504a238c70SJohn Marino int k, l;
514a238c70SJohn Marino
524a238c70SJohn Marino MPFR_ASSERTN ((size_t) m < sizeof (long) * CHAR_BIT - 1);
534a238c70SJohn Marino
544a238c70SJohn Marino S = Q + (m+1);
554a238c70SJohn Marino ptoj = Q + 2*(m+1); /* ptoj[i] = mantissa^(2^i) */
564a238c70SJohn Marino log2_nb_terms = mult + (m+1);
574a238c70SJohn Marino
584a238c70SJohn Marino /* Normalize p */
594a238c70SJohn Marino MPFR_ASSERTD (mpz_cmp_ui (p, 0) != 0);
604a238c70SJohn Marino n = mpz_scan1 (p, 0); /* number of trailing zeros in p */
614a238c70SJohn Marino mpz_tdiv_q_2exp (p, p, n);
624a238c70SJohn Marino r -= n; /* since |p/2^r| < 1 and p >= 1, r >= 1 */
634a238c70SJohn Marino
644a238c70SJohn Marino /* Set initial var */
654a238c70SJohn Marino mpz_set (ptoj[0], p);
664a238c70SJohn Marino for (k = 1; k < m; k++)
674a238c70SJohn Marino mpz_mul (ptoj[k], ptoj[k-1], ptoj[k-1]); /* ptoj[k] = p^(2^k) */
684a238c70SJohn Marino mpz_set_ui (Q[0], 1);
694a238c70SJohn Marino mpz_set_ui (S[0], 1);
704a238c70SJohn Marino k = 0;
714a238c70SJohn Marino mult[0] = 0; /* the multiplier P[k]/Q[k] for the remaining terms
724a238c70SJohn Marino satisfies P[k]/Q[k] <= 2^(-mult[k]) */
734a238c70SJohn Marino log2_nb_terms[0] = 0; /* log2(#terms) [exact in 1st loop where 2^k] */
744a238c70SJohn Marino prec_i_have = 0;
754a238c70SJohn Marino
764a238c70SJohn Marino /* Main Loop */
774a238c70SJohn Marino n = 1UL << m;
784a238c70SJohn Marino for (i = 1; (prec_i_have < precy) && (i < n); i++)
794a238c70SJohn Marino {
804a238c70SJohn Marino /* invariant: Q[0]*Q[1]*...*Q[k] equals i! */
814a238c70SJohn Marino k++;
824a238c70SJohn Marino log2_nb_terms[k] = 0; /* 1 term */
834a238c70SJohn Marino mpz_set_ui (Q[k], i + 1);
844a238c70SJohn Marino mpz_set_ui (S[k], i + 1);
854a238c70SJohn Marino j = i + 1; /* we have computed j = i+1 terms so far */
864a238c70SJohn Marino l = 0;
874a238c70SJohn Marino while ((j & 1) == 0) /* combine and reduce */
884a238c70SJohn Marino {
894a238c70SJohn Marino /* invariant: S[k] corresponds to 2^l consecutive terms */
904a238c70SJohn Marino mpz_mul (S[k], S[k], ptoj[l]);
914a238c70SJohn Marino mpz_mul (S[k-1], S[k-1], Q[k]);
924a238c70SJohn Marino /* Q[k] corresponds to 2^l consecutive terms too.
934a238c70SJohn Marino Since it does not contains the factor 2^(r*2^l),
944a238c70SJohn Marino when going from l to l+1 we need to multiply
954a238c70SJohn Marino by 2^(r*2^(l+1))/2^(r*2^l) = 2^(r*2^l) */
964a238c70SJohn Marino mpz_mul_2exp (S[k-1], S[k-1], r << l);
974a238c70SJohn Marino mpz_add (S[k-1], S[k-1], S[k]);
984a238c70SJohn Marino mpz_mul (Q[k-1], Q[k-1], Q[k]);
994a238c70SJohn Marino log2_nb_terms[k-1] ++; /* number of terms in S[k-1]
1004a238c70SJohn Marino is a power of 2 by construction */
1014a238c70SJohn Marino MPFR_MPZ_SIZEINBASE2 (prec_i_have, Q[k]);
1024a238c70SJohn Marino MPFR_MPZ_SIZEINBASE2 (prec_ptoj, ptoj[l]);
1034a238c70SJohn Marino mult[k-1] += prec_i_have + (r << l) - prec_ptoj - 1;
1044a238c70SJohn Marino prec_i_have = mult[k] = mult[k-1];
1054a238c70SJohn Marino /* since mult[k] >= mult[k-1] + nbits(Q[k]),
1064a238c70SJohn Marino we have Q[0]*...*Q[k] <= 2^mult[k] = 2^prec_i_have */
1074a238c70SJohn Marino l ++;
1084a238c70SJohn Marino j >>= 1;
1094a238c70SJohn Marino k --;
1104a238c70SJohn Marino }
1114a238c70SJohn Marino }
1124a238c70SJohn Marino
1134a238c70SJohn Marino /* accumulate all products in S[0] and Q[0]. Warning: contrary to above,
1144a238c70SJohn Marino here we do not have log2_nb_terms[k-1] = log2_nb_terms[k]+1. */
1154a238c70SJohn Marino l = 0; /* number of accumulated terms in the right part S[k]/Q[k] */
1164a238c70SJohn Marino while (k > 0)
1174a238c70SJohn Marino {
1184a238c70SJohn Marino j = log2_nb_terms[k-1];
1194a238c70SJohn Marino mpz_mul (S[k], S[k], ptoj[j]);
1204a238c70SJohn Marino mpz_mul (S[k-1], S[k-1], Q[k]);
1214a238c70SJohn Marino l += 1 << log2_nb_terms[k];
1224a238c70SJohn Marino mpz_mul_2exp (S[k-1], S[k-1], r * l);
1234a238c70SJohn Marino mpz_add (S[k-1], S[k-1], S[k]);
1244a238c70SJohn Marino mpz_mul (Q[k-1], Q[k-1], Q[k]);
1254a238c70SJohn Marino k--;
1264a238c70SJohn Marino }
1274a238c70SJohn Marino
1284a238c70SJohn Marino /* Q[0] now equals i! */
1294a238c70SJohn Marino MPFR_MPZ_SIZEINBASE2 (prec_i_have, S[0]);
1304a238c70SJohn Marino diff = (mpfr_exp_t) prec_i_have - 2 * (mpfr_exp_t) precy;
1314a238c70SJohn Marino expo = diff;
1324a238c70SJohn Marino if (diff >= 0)
1334a238c70SJohn Marino mpz_fdiv_q_2exp (S[0], S[0], diff);
1344a238c70SJohn Marino else
1354a238c70SJohn Marino mpz_mul_2exp (S[0], S[0], -diff);
1364a238c70SJohn Marino
1374a238c70SJohn Marino MPFR_MPZ_SIZEINBASE2 (prec_i_have, Q[0]);
1384a238c70SJohn Marino diff = (mpfr_exp_t) prec_i_have - (mpfr_prec_t) precy;
1394a238c70SJohn Marino expo -= diff;
1404a238c70SJohn Marino if (diff > 0)
1414a238c70SJohn Marino mpz_fdiv_q_2exp (Q[0], Q[0], diff);
1424a238c70SJohn Marino else
1434a238c70SJohn Marino mpz_mul_2exp (Q[0], Q[0], -diff);
1444a238c70SJohn Marino
1454a238c70SJohn Marino mpz_tdiv_q (S[0], S[0], Q[0]);
1464a238c70SJohn Marino mpfr_set_z (y, S[0], MPFR_RNDD);
1474a238c70SJohn Marino MPFR_SET_EXP (y, MPFR_GET_EXP (y) + expo - r * (i - 1) );
1484a238c70SJohn Marino }
1494a238c70SJohn Marino
1504a238c70SJohn Marino #define shift (GMP_NUMB_BITS/2)
1514a238c70SJohn Marino
1524a238c70SJohn Marino int
mpfr_exp_3(mpfr_ptr y,mpfr_srcptr x,mpfr_rnd_t rnd_mode)1534a238c70SJohn Marino mpfr_exp_3 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
1544a238c70SJohn Marino {
1554a238c70SJohn Marino mpfr_t t, x_copy, tmp;
1564a238c70SJohn Marino mpz_t uk;
1574a238c70SJohn Marino mpfr_exp_t ttt, shift_x;
1584a238c70SJohn Marino unsigned long twopoweri;
1594a238c70SJohn Marino mpz_t *P;
1604a238c70SJohn Marino mpfr_prec_t *mult;
1614a238c70SJohn Marino int i, k, loop;
1624a238c70SJohn Marino int prec_x;
1634a238c70SJohn Marino mpfr_prec_t realprec, Prec;
1644a238c70SJohn Marino int iter;
1654a238c70SJohn Marino int inexact = 0;
1664a238c70SJohn Marino MPFR_SAVE_EXPO_DECL (expo);
1674a238c70SJohn Marino MPFR_ZIV_DECL (ziv_loop);
1684a238c70SJohn Marino
1694a238c70SJohn Marino MPFR_LOG_FUNC
1704a238c70SJohn Marino (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec(x), mpfr_log_prec, x, rnd_mode),
1714a238c70SJohn Marino ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec(y), mpfr_log_prec, y,
1724a238c70SJohn Marino inexact));
1734a238c70SJohn Marino
1744a238c70SJohn Marino MPFR_SAVE_EXPO_MARK (expo);
1754a238c70SJohn Marino
1764a238c70SJohn Marino /* decompose x */
1774a238c70SJohn Marino /* we first write x = 1.xxxxxxxxxxxxx
1784a238c70SJohn Marino ----- k bits -- */
1794a238c70SJohn Marino prec_x = MPFR_INT_CEIL_LOG2 (MPFR_PREC (x)) - MPFR_LOG2_GMP_NUMB_BITS;
1804a238c70SJohn Marino if (prec_x < 0)
1814a238c70SJohn Marino prec_x = 0;
1824a238c70SJohn Marino
1834a238c70SJohn Marino ttt = MPFR_GET_EXP (x);
1844a238c70SJohn Marino mpfr_init2 (x_copy, MPFR_PREC(x));
1854a238c70SJohn Marino mpfr_set (x_copy, x, MPFR_RNDD);
1864a238c70SJohn Marino
1874a238c70SJohn Marino /* we shift to get a number less than 1 */
1884a238c70SJohn Marino if (ttt > 0)
1894a238c70SJohn Marino {
1904a238c70SJohn Marino shift_x = ttt;
1914a238c70SJohn Marino mpfr_div_2ui (x_copy, x, ttt, MPFR_RNDN);
1924a238c70SJohn Marino ttt = MPFR_GET_EXP (x_copy);
1934a238c70SJohn Marino }
1944a238c70SJohn Marino else
1954a238c70SJohn Marino shift_x = 0;
1964a238c70SJohn Marino MPFR_ASSERTD (ttt <= 0);
1974a238c70SJohn Marino
1984a238c70SJohn Marino /* Init prec and vars */
1994a238c70SJohn Marino realprec = MPFR_PREC (y) + MPFR_INT_CEIL_LOG2 (prec_x + MPFR_PREC (y));
2004a238c70SJohn Marino Prec = realprec + shift + 2 + shift_x;
2014a238c70SJohn Marino mpfr_init2 (t, Prec);
2024a238c70SJohn Marino mpfr_init2 (tmp, Prec);
2034a238c70SJohn Marino mpz_init (uk);
2044a238c70SJohn Marino
2054a238c70SJohn Marino /* Main loop */
2064a238c70SJohn Marino MPFR_ZIV_INIT (ziv_loop, realprec);
2074a238c70SJohn Marino for (;;)
2084a238c70SJohn Marino {
2094a238c70SJohn Marino int scaled = 0;
2104a238c70SJohn Marino MPFR_BLOCK_DECL (flags);
2114a238c70SJohn Marino
2124a238c70SJohn Marino k = MPFR_INT_CEIL_LOG2 (Prec) - MPFR_LOG2_GMP_NUMB_BITS;
2134a238c70SJohn Marino
2144a238c70SJohn Marino /* now we have to extract */
2154a238c70SJohn Marino twopoweri = GMP_NUMB_BITS;
2164a238c70SJohn Marino
2174a238c70SJohn Marino /* Allocate tables */
2184a238c70SJohn Marino P = (mpz_t*) (*__gmp_allocate_func) (3*(k+2)*sizeof(mpz_t));
2194a238c70SJohn Marino for (i = 0; i < 3*(k+2); i++)
2204a238c70SJohn Marino mpz_init (P[i]);
2214a238c70SJohn Marino mult = (mpfr_prec_t*) (*__gmp_allocate_func) (2*(k+2)*sizeof(mpfr_prec_t));
2224a238c70SJohn Marino
2234a238c70SJohn Marino /* Particular case for i==0 */
2244a238c70SJohn Marino mpfr_extract (uk, x_copy, 0);
2254a238c70SJohn Marino MPFR_ASSERTD (mpz_cmp_ui (uk, 0) != 0);
2264a238c70SJohn Marino mpfr_exp_rational (tmp, uk, shift + twopoweri - ttt, k + 1, P, mult);
2274a238c70SJohn Marino for (loop = 0; loop < shift; loop++)
2284a238c70SJohn Marino mpfr_sqr (tmp, tmp, MPFR_RNDD);
2294a238c70SJohn Marino twopoweri *= 2;
2304a238c70SJohn Marino
2314a238c70SJohn Marino /* General case */
2324a238c70SJohn Marino iter = (k <= prec_x) ? k : prec_x;
2334a238c70SJohn Marino for (i = 1; i <= iter; i++)
2344a238c70SJohn Marino {
2354a238c70SJohn Marino mpfr_extract (uk, x_copy, i);
2364a238c70SJohn Marino if (MPFR_LIKELY (mpz_cmp_ui (uk, 0) != 0))
2374a238c70SJohn Marino {
2384a238c70SJohn Marino mpfr_exp_rational (t, uk, twopoweri - ttt, k - i + 1, P, mult);
2394a238c70SJohn Marino mpfr_mul (tmp, tmp, t, MPFR_RNDD);
2404a238c70SJohn Marino }
2414a238c70SJohn Marino MPFR_ASSERTN (twopoweri <= LONG_MAX/2);
2424a238c70SJohn Marino twopoweri *=2;
2434a238c70SJohn Marino }
2444a238c70SJohn Marino
2454a238c70SJohn Marino /* Clear tables */
2464a238c70SJohn Marino for (i = 0; i < 3*(k+2); i++)
2474a238c70SJohn Marino mpz_clear (P[i]);
2484a238c70SJohn Marino (*__gmp_free_func) (P, 3*(k+2)*sizeof(mpz_t));
2494a238c70SJohn Marino (*__gmp_free_func) (mult, 2*(k+2)*sizeof(mpfr_prec_t));
2504a238c70SJohn Marino
2514a238c70SJohn Marino if (shift_x > 0)
2524a238c70SJohn Marino {
2534a238c70SJohn Marino MPFR_BLOCK (flags, {
2544a238c70SJohn Marino for (loop = 0; loop < shift_x - 1; loop++)
2554a238c70SJohn Marino mpfr_sqr (tmp, tmp, MPFR_RNDD);
2564a238c70SJohn Marino mpfr_sqr (t, tmp, MPFR_RNDD);
2574a238c70SJohn Marino } );
2584a238c70SJohn Marino
2594a238c70SJohn Marino if (MPFR_UNLIKELY (MPFR_OVERFLOW (flags)))
2604a238c70SJohn Marino {
2614a238c70SJohn Marino /* tmp <= exact result, so that it is a real overflow. */
2624a238c70SJohn Marino inexact = mpfr_overflow (y, rnd_mode, 1);
2634a238c70SJohn Marino MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW);
2644a238c70SJohn Marino break;
2654a238c70SJohn Marino }
2664a238c70SJohn Marino
2674a238c70SJohn Marino if (MPFR_UNLIKELY (MPFR_UNDERFLOW (flags)))
2684a238c70SJohn Marino {
2694a238c70SJohn Marino /* This may be a spurious underflow. So, let's scale
2704a238c70SJohn Marino the result. */
2714a238c70SJohn Marino mpfr_mul_2ui (tmp, tmp, 1, MPFR_RNDD); /* no overflow, exact */
2724a238c70SJohn Marino mpfr_sqr (t, tmp, MPFR_RNDD);
2734a238c70SJohn Marino if (MPFR_IS_ZERO (t))
2744a238c70SJohn Marino {
2754a238c70SJohn Marino /* approximate result < 2^(emin - 3), thus
2764a238c70SJohn Marino exact result < 2^(emin - 2). */
2774a238c70SJohn Marino inexact = mpfr_underflow (y, (rnd_mode == MPFR_RNDN) ?
2784a238c70SJohn Marino MPFR_RNDZ : rnd_mode, 1);
2794a238c70SJohn Marino MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_UNDERFLOW);
2804a238c70SJohn Marino break;
2814a238c70SJohn Marino }
2824a238c70SJohn Marino scaled = 1;
2834a238c70SJohn Marino }
2844a238c70SJohn Marino }
2854a238c70SJohn Marino
2864a238c70SJohn Marino if (mpfr_can_round (shift_x > 0 ? t : tmp, realprec, MPFR_RNDD, MPFR_RNDZ,
2874a238c70SJohn Marino MPFR_PREC(y) + (rnd_mode == MPFR_RNDN)))
2884a238c70SJohn Marino {
2894a238c70SJohn Marino inexact = mpfr_set (y, shift_x > 0 ? t : tmp, rnd_mode);
2904a238c70SJohn Marino if (MPFR_UNLIKELY (scaled && MPFR_IS_PURE_FP (y)))
2914a238c70SJohn Marino {
2924a238c70SJohn Marino int inex2;
2934a238c70SJohn Marino mpfr_exp_t ey;
2944a238c70SJohn Marino
2954a238c70SJohn Marino /* The result has been scaled and needs to be corrected. */
2964a238c70SJohn Marino ey = MPFR_GET_EXP (y);
2974a238c70SJohn Marino inex2 = mpfr_mul_2si (y, y, -2, rnd_mode);
2984a238c70SJohn Marino if (inex2) /* underflow */
2994a238c70SJohn Marino {
3004a238c70SJohn Marino if (rnd_mode == MPFR_RNDN && inexact < 0 &&
3014a238c70SJohn Marino MPFR_IS_ZERO (y) && ey == __gmpfr_emin + 1)
3024a238c70SJohn Marino {
3034a238c70SJohn Marino /* Double rounding case: in MPFR_RNDN, the scaled
3044a238c70SJohn Marino result has been rounded downward to 2^emin.
3054a238c70SJohn Marino As the exact result is > 2^(emin - 2), correct
3064a238c70SJohn Marino rounding must be done upward. */
3074a238c70SJohn Marino /* TODO: make sure in coverage tests that this line
3084a238c70SJohn Marino is reached. */
3094a238c70SJohn Marino inexact = mpfr_underflow (y, MPFR_RNDU, 1);
3104a238c70SJohn Marino }
3114a238c70SJohn Marino else
3124a238c70SJohn Marino {
3134a238c70SJohn Marino /* No double rounding. */
3144a238c70SJohn Marino inexact = inex2;
3154a238c70SJohn Marino }
3164a238c70SJohn Marino MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_UNDERFLOW);
3174a238c70SJohn Marino }
3184a238c70SJohn Marino }
3194a238c70SJohn Marino break;
3204a238c70SJohn Marino }
3214a238c70SJohn Marino
3224a238c70SJohn Marino MPFR_ZIV_NEXT (ziv_loop, realprec);
3234a238c70SJohn Marino Prec = realprec + shift + 2 + shift_x;
3244a238c70SJohn Marino mpfr_set_prec (t, Prec);
3254a238c70SJohn Marino mpfr_set_prec (tmp, Prec);
3264a238c70SJohn Marino }
3274a238c70SJohn Marino MPFR_ZIV_FREE (ziv_loop);
3284a238c70SJohn Marino
3294a238c70SJohn Marino mpz_clear (uk);
3304a238c70SJohn Marino mpfr_clear (tmp);
3314a238c70SJohn Marino mpfr_clear (t);
3324a238c70SJohn Marino mpfr_clear (x_copy);
3334a238c70SJohn Marino MPFR_SAVE_EXPO_FREE (expo);
3344a238c70SJohn Marino return inexact;
3354a238c70SJohn Marino }
336