1
2 /*
3 * M_APM - mapmpwr2.c
4 *
5 * Copyright (C) 2002 - 2007 Michael C. Ring
6 *
7 * Permission to use, copy, and distribute this software and its
8 * documentation for any purpose with or without fee is hereby granted,
9 * provided that the above copyright notice appear in all copies and
10 * that both that copyright notice and this permission notice appear
11 * in supporting documentation.
12 *
13 * Permission to modify the software is granted. Permission to distribute
14 * the modified code is granted. Modifications are to be distributed by
15 * using the file 'license.txt' as a template to modify the file header.
16 * 'license.txt' is available in the official MAPM distribution.
17 *
18 * This software is provided "as is" without express or implied warranty.
19 */
20
21 /*
22 * $Id: mapmpwr2.c,v 1.4 2007/12/03 01:56:21 mike Exp $
23 *
24 * This file contains the Integer Power function and the result
25 * is NOT ROUNDED. The exponent must be an integer >= zero.
26 *
27 * This will typically be used in an application where full integer
28 * precision is required to be maintained.
29 *
30 * $Log: mapmpwr2.c,v $
31 * Revision 1.4 2007/12/03 01:56:21 mike
32 * Update license
33 *
34 * Revision 1.3 2003/07/21 20:38:06 mike
35 * Modify error messages to be in a consistent format.
36 *
37 * Revision 1.2 2003/03/31 21:51:23 mike
38 * call generic error handling function
39 *
40 * Revision 1.1 2002/11/03 21:02:04 mike
41 * Initial revision
42 */
43
44 #include "m_apm_lc.h"
45
46 /****************************************************************************/
m_apm_integer_pow_nr(M_APM rr,M_APM aa,int mexp)47 void m_apm_integer_pow_nr(M_APM rr, M_APM aa, int mexp)
48 {
49 M_APM tmp0, tmpy, tmpz;
50 int nexp, ii;
51
52 if (mexp == 0)
53 {
54 m_apm_copy(rr, MM_One);
55 return;
56 }
57 else
58 {
59 if (mexp < 0)
60 {
61 M_apm_log_error_msg(M_APM_RETURN,
62 "\'m_apm_integer_pow_nr\', Negative exponent");
63
64 M_set_to_zero(rr);
65 return;
66 }
67 }
68
69 if (mexp == 1)
70 {
71 m_apm_copy(rr, aa);
72 return;
73 }
74
75 if (mexp == 2)
76 {
77 m_apm_multiply(rr, aa, aa);
78 return;
79 }
80
81 nexp = mexp;
82
83 if (aa->m_apm_sign == 0)
84 {
85 M_set_to_zero(rr);
86 return;
87 }
88
89 tmp0 = M_get_stack_var();
90 tmpy = M_get_stack_var();
91 tmpz = M_get_stack_var();
92
93 m_apm_copy(tmpy, MM_One);
94 m_apm_copy(tmpz, aa);
95
96 while (TRUE)
97 {
98 ii = nexp & 1;
99 nexp = nexp >> 1;
100
101 if (ii != 0) /* exponent -was- odd */
102 {
103 m_apm_multiply(tmp0, tmpy, tmpz);
104
105 if (nexp == 0)
106 break;
107
108 m_apm_copy(tmpy, tmp0);
109 }
110
111 m_apm_multiply(tmp0, tmpz, tmpz);
112 m_apm_copy(tmpz, tmp0);
113 }
114
115 m_apm_copy(rr, tmp0);
116
117 M_restore_stack(3);
118 }
119 /****************************************************************************/
120
121