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 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 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