1 2 /* 3 * M_APM - mapmipwr.c 4 * 5 * Copyright (C) 1999 - 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: mapmipwr.c,v 1.6 2007/12/03 01:54:39 mike Exp $ 23 * 24 * This file contains the Integer Power function. 25 * 26 * $Log: mapmipwr.c,v $ 27 * Revision 1.6 2007/12/03 01:54:39 mike 28 * Update license 29 * 30 * Revision 1.5 2002/11/03 21:10:32 mike 31 * Updated function parameters to use the modern style 32 * 33 * Revision 1.4 2000/09/23 19:46:04 mike 34 * change divide call to reciprocal 35 * 36 * Revision 1.3 2000/05/24 17:03:35 mike 37 * return 1 when input is 0^0 38 * 39 * Revision 1.2 1999/09/18 01:34:35 mike 40 * minor tweaks 41 * 42 * Revision 1.1 1999/09/18 01:33:09 mike 43 * Initial revision 44 */ 45 46 #include "m_apm_lc.h" 47 48 /****************************************************************************/ 49 void m_apm_integer_pow(M_APM rr, int places, M_APM aa, int mexp) 50 { 51 M_APM tmp0, tmpy, tmpz; 52 int nexp, ii, signflag, local_precision; 53 54 if (mexp == 0) 55 { 56 m_apm_copy(rr, MM_One); 57 return; 58 } 59 else 60 { 61 if (mexp > 0) 62 { 63 signflag = 0; 64 nexp = mexp; 65 } 66 else 67 { 68 signflag = 1; 69 nexp = -mexp; 70 } 71 } 72 73 if (aa->m_apm_sign == 0) 74 { 75 M_set_to_zero(rr); 76 return; 77 } 78 79 tmp0 = M_get_stack_var(); 80 tmpy = M_get_stack_var(); 81 tmpz = M_get_stack_var(); 82 83 local_precision = places + 8; 84 85 m_apm_copy(tmpy, MM_One); 86 m_apm_copy(tmpz, aa); 87 88 while (TRUE) 89 { 90 ii = nexp & 1; 91 nexp = nexp >> 1; 92 93 if (ii != 0) /* exponent -was- odd */ 94 { 95 m_apm_multiply(tmp0, tmpy, tmpz); 96 m_apm_round(tmpy, local_precision, tmp0); 97 98 if (nexp == 0) 99 break; 100 } 101 102 m_apm_multiply(tmp0, tmpz, tmpz); 103 m_apm_round(tmpz, local_precision, tmp0); 104 } 105 106 if (signflag) 107 { 108 m_apm_reciprocal(rr, places, tmpy); 109 } 110 else 111 { 112 m_apm_round(rr, places, tmpy); 113 } 114 115 M_restore_stack(3); 116 } 117 /****************************************************************************/ 118