1 /* $OpenBSD: epow.c,v 1.1 2011/07/02 18:11:01 martynas Exp $ */
2
3 /*
4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5 *
6 * Permission to use, copy, modify, and distribute this software for any
7 * purpose with or without fee is hereby granted, provided that the above
8 * copyright notice and this permission notice appear in all copies.
9 *
10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17 */
18
19 /* epow.c */
20 /* power function: z = x**y */
21 /* by Stephen L. Moshier. */
22
23
24 #include "ehead.h"
25
26 extern int rndprc;
27 void epowi();
28
epow(x,y,z)29 void epow( x, y, z )
30 unsigned short *x, *y, *z;
31 {
32 unsigned short w[NE];
33 int rndsav;
34 long li;
35
36 efloor( y, w );
37 if( ecmp(y,w) == 0 )
38 {
39 eifrac( y, &li, w );
40 if( li < 0 )
41 li = -li;
42 if( (li < 0x7fffffff) && (li != 0x80000000) )
43 {
44 epowi( x, y, z );
45 return;
46 }
47 }
48 /* z = exp( y * log(x) ) */
49 rndsav = rndprc;
50 rndprc = NBITS;
51 elog( x, w );
52 emul( y, w, w );
53 eexp( w, z );
54 rndprc = rndsav;
55 emul( eone, z, z );
56 }
57
58
59 /* y is integer valued. */
60
epowi(x,y,z)61 void epowi( x, y, z )
62 unsigned short x[], y[], z[];
63 {
64 unsigned short w[NE];
65 long li, lx;
66 unsigned long lu;
67 int rndsav;
68 unsigned short signx;
69 /* unsigned short signy; */
70
71 rndsav = rndprc;
72 eifrac( y, &li, w );
73 if( li < 0 )
74 lx = -li;
75 else
76 lx = li;
77
78 if( (lx == 0x7fffffff) || (lx == 0x80000000) )
79 {
80 epow( x, y, z );
81 goto done;
82 }
83
84 if( (x[NE-1] & (unsigned short )0x7fff) == 0 )
85 {
86 if( li == 0 )
87 {
88 emov( eone, z );
89 return;
90 }
91 else if( li < 0 )
92 {
93 einfin( z );
94 return;
95 }
96 else
97 {
98 eclear( z );
99 return;
100 }
101 }
102
103 if( li == 0L )
104 {
105 emov( eone, z );
106 return;
107 }
108
109 emov( x, w );
110 signx = w[NE-1] & (unsigned short )0x8000;
111 w[NE-1] &= (unsigned short )0x7fff;
112
113 /* Overflow detection */
114 /*
115 lx = li * (w[NE-1] - 0x3fff);
116 if( lx > 16385L )
117 {
118 einfin( z );
119 mtherr( "epowi", OVERFLOW );
120 goto done;
121 }
122 if( lx < -16450L )
123 {
124 eclear( z );
125 return;
126 }
127 */
128 rndprc = NBITS;
129
130 if( li < 0 )
131 {
132 lu = (unsigned int )( -li );
133 /* signy = 0xffff;*/
134 ediv( w, eone, w );
135 }
136 else
137 {
138 lu = (unsigned int )li;
139 /* signy = 0;*/
140 }
141
142 /* First bit of the power */
143 if( lu & 1 )
144 {
145 emov( w, z );
146 }
147 else
148 {
149 emov( eone, z );
150 signx = 0;
151 }
152
153
154 lu >>= 1;
155 while( lu != 0L )
156 {
157 emul( w, w, w ); /* arg to the 2-to-the-kth power */
158 if( lu & 1L ) /* if that bit is set, then include in product */
159 emul( w, z, z );
160 lu >>= 1;
161 }
162
163
164 done:
165
166 if( signx )
167 eneg( z ); /* odd power of negative number */
168
169 /*
170 if( signy )
171 {
172 if( ecmp( z, ezero ) != 0 )
173 {
174 ediv( z, eone, z );
175 }
176 else
177 {
178 einfin( z );
179 printf( "epowi OVERFLOW\n" );
180 }
181 }
182 */
183 rndprc = rndsav;
184 emul( eone, z, z );
185 }
186
187
188