qpow.cpp

   1 // Rational power function
   2
   3 #include "stdafx.h"
   4 #include "defs.h"
   5
   6 static void qpowf(void);
   7 static void normalize_angle(void);
   8 static int is_small_integer(U *);
   9
  10 void
  11 qpow()
  12 {
  13         save();
  14         qpowf();
  15         restore();
  16 }
  17
  18 #define BASE p1
  19 #define EXPO p2
  20
  21 static void
  22 qpowf(void)
  23 {
  24         int expo;
  25         unsigned int a, b, *t, *x, *y;
  26
  27         EXPO = pop();
  28         BASE = pop();
  29
  30         // if base is 1 or exponent is 0 then return 1
  31
  32         if (isplusone(BASE) || iszero(EXPO)) {
  33                 push_integer(1);
  34                 return;
  35         }
  36
  37         // if base is zero then return 0
  38
  39         if (iszero(BASE)) {
  40                 if (isnegativenumber(EXPO))
  41                         stop("divide by zero");
  42                 push(zero);
  43                 return;
  44         }
  45
  46         // if exponent is 1 then return base
  47
  48         if (isplusone(EXPO)) {
  49                 push(BASE);
  50                 return;
  51         }
  52
  53         // if exponent is integer then power
  54
  55         if (isinteger(EXPO)) {
  56                 push(EXPO);
  57                 expo = pop_integer();
  58                 if (expo == (int) 0x80000000) {
  59                         // expo greater than 32 bits
  60                         push_symbol(POWER);
  61                         push(BASE);
  62                         push(EXPO);
  63                         list(3);
  64                         return;
  65                 }
  66                 x = mpow(BASE->u.q.a, abs(expo));
  67                 y = mpow(BASE->u.q.b, abs(expo));
  68                 if (expo < 0) {
  69                         t = x;
  70                         x = y;
  71                         y = t;
  72                         MSIGN(x) = MSIGN(y);
  73                         MSIGN(y) = 1;
  74                 }
  75                 p3 = alloc();
  76                 p3->k = NUM;
  77                 p3->u.q.a = x;
  78                 p3->u.q.b = y;
  79                 push(p3);
  80                 return;
  81         }
  82
  83         // from here on out the exponent is NOT an integer
  84
  85         // if base is -1 then normalize polar angle
  86
  87         if (isminusone(BASE)) {
  88                 push(EXPO);
  89                 normalize_angle();
  90                 return;
  91         }
  92
  93         // if base is negative then (-N)^M -> N^M * (-1)^M
  94
  95         if (isnegativenumber(BASE)) {
  96                 push(BASE);
  97                 negate();
  98                 push(EXPO);
  99                 qpow();
 100                 push_integer(-1);
 101                 push(EXPO);
 102                 qpow();
 103                 multiply();
 104                 return;
 105         }
 106
 107         // if BASE is not an integer then power numerator and denominator
 108
 109         if (!isinteger(BASE)) {
 110                 push(BASE);
 111                 mp_numerator();
 112                 push(EXPO);
 113                 qpow();
 114                 push(BASE);
 115                 mp_denominator();
 116                 push(EXPO);
 117                 negate();
 118                 qpow();
 119                 multiply();
 120                 return;
 121         }
 122
 123         // At this point BASE is a positive integer.
 124
 125         // If BASE is small then factor it.
 126
 127         if (is_small_integer(BASE)) {
 128                 push(BASE);
 129                 push(EXPO);
 130                 quickfactor();
 131                 return;
 132         }
 133
 134         // At this point BASE is a positive integer and EXPO is not an integer.
 135
 136         if (MLENGTH(EXPO->u.q.a) > 1 || MLENGTH(EXPO->u.q.b) > 1) {
 137                 push_symbol(POWER);
 138                 push(BASE);
 139                 push(EXPO);
 140                 list(3);
 141                 return;
 142         }
 143
 144         a = EXPO->u.q.a[0];
 145         b = EXPO->u.q.b[0];
 146
 147         x = mroot(BASE->u.q.a, b);
 148
 149         if (x == 0) {
 150                 push_symbol(POWER);
 151                 push(BASE);
 152                 push(EXPO);
 153                 list(3);
 154                 return;
 155         }
 156
 157         y = mpow(x, a);
 158
 159         mfree(x);
 160
 161         p3 = alloc();
 162
 163         p3->k = NUM;
 164
 165         if (MSIGN(EXPO->u.q.a) == -1) {
 166                 p3->u.q.a = mint(1);
 167                 p3->u.q.b = y;
 168         } else {
 169                 p3->u.q.a = y;
 170                 p3->u.q.b = mint(1);
 171         }
 172
 173         push(p3);
 174 }
 175
 176 //-----------------------------------------------------------------------------
 177 //
 178 //      Normalize the angle of unit imaginary, i.e. (-1) ^ N
 179 //
 180 //      Input:          N on stack (must be rational, not float)
 181 //
 182 //      Output:         Result on stack
 183 //
 184 //      Note:
 185 //
 186 //      n = q * d + r
 187 //
 188 //      Example:
 189 //                                              n       d       q       r
 190 //
 191 //      (-1)^(8/3)      ->       (-1)^(2/3)     8       3       2       2
 192 //      (-1)^(7/3)      ->       (-1)^(1/3)     7       3       2       1
 193 //      (-1)^(5/3)      ->      -(-1)^(2/3)     5       3       1       2
 194 //      (-1)^(4/3)      ->      -(-1)^(1/3)     4       3       1       1
 195 //      (-1)^(2/3)      ->       (-1)^(2/3)     2       3       0       2
 196 //      (-1)^(1/3)      ->       (-1)^(1/3)     1       3       0       1
 197 //
 198 //      (-1)^(-1/3)     ->      -(-1)^(2/3)     -1      3       -1      2
 199 //      (-1)^(-2/3)     ->      -(-1)^(1/3)     -2      3       -1      1
 200 //      (-1)^(-4/3)     ->       (-1)^(2/3)     -4      3       -2      2
 201 //      (-1)^(-5/3)     ->       (-1)^(1/3)     -5      3       -2      1
 202 //      (-1)^(-7/3)     ->      -(-1)^(2/3)     -7      3       -3      2
 203 //      (-1)^(-8/3)     ->      -(-1)^(1/3)     -8      3       -3      1
 204 //
 205 //-----------------------------------------------------------------------------
 206
 207 #define A p1
 208 #define Q p2
 209 #define R p3
 210
 211 static void
 212 normalize_angle(void)
 213 {
 214         save();
 215
 216         A = pop();
 217
 218         // integer exponent?
 219
 220         if (isinteger(A)) {
 221                 if (A->u.q.a[0] & 1)
 222                         push_integer(-1); // odd exponent
 223                 else
 224                         push_integer(1); // even exponent
 225                 restore();
 226                 return;
 227         }
 228
 229         // floor
 230
 231         push(A);
 232         bignum_truncate();
 233         Q = pop();
 234
 235         if (isnegativenumber(A)) {
 236                 push(Q);
 237                 push_integer(-1);
 238                 add();
 239                 Q = pop();
 240         }
 241
 242         // remainder (always positive)
 243
 244         push(A);
 245         push(Q);
 246         subtract();
 247         R = pop();
 248
 249         // remainder becomes new angle
 250
 251         push_symbol(POWER);
 252         push_integer(-1);
 253         push(R);
 254         list(3);
 255
 256         // negate if quotient is odd
 257
 258         if (Q->u.q.a[0] & 1)
 259                 negate();
 260
 261         restore();
 262 }
 263
 264 static int
 265 is_small_integer(U *p)
 266 {
 267         if (isinteger(p) && MLENGTH(p->u.q.a) == 1 && (p->u.q.a[0] & 0x80000000) == 0)
 268                 return 1;
 269         else
 270                 return 0;
 271 }