// Rational power function #include "stdafx.h" #include "defs.h" static void qpowf(void); static void normalize_angle(void); static int is_small_integer(U *); void qpow() { save(); qpowf(); restore(); } #define BASE p1 #define EXPO p2 static void qpowf(void) { int expo; unsigned int a, b, *t, *x, *y; EXPO = pop(); BASE = pop(); // if base is 1 or exponent is 0 then return 1 if (isplusone(BASE) || iszero(EXPO)) { push_integer(1); return; } // if base is zero then return 0 if (iszero(BASE)) { if (isnegativenumber(EXPO)) stop("divide by zero"); push(_zero); return; } // if exponent is 1 then return base if (isplusone(EXPO)) { push(BASE); return; } // if exponent is integer then power if (isinteger(EXPO)) { push(EXPO); expo = pop_integer(); if (expo == (int) 0x80000000) { // expo greater than 32 bits push_symbol(POWER); push(BASE); push(EXPO); list(3); return; } x = mpow(BASE->u.q.a, abs(expo)); y = mpow(BASE->u.q.b, abs(expo)); if (expo < 0) { t = x; x = y; y = t; MSIGN(x) = MSIGN(y); MSIGN(y) = 1; } p3 = alloc(); p3->k = NUM; p3->u.q.a = x; p3->u.q.b = y; push(p3); return; } // from here on out the exponent is NOT an integer // if base is -1 then normalize polar angle if (isminusone(BASE)) { push(EXPO); if (conjugating) negate(); normalize_angle(); return; } // if base is negative then (-N)^M -> N^M * (-1)^M if (isnegativenumber(BASE)) { push(BASE); negate(); push(EXPO); power_numbers(); push_integer(-1); push(EXPO); power_numbers(); multiply(); return; } // if BASE is not an integer then power numerator and denominator if (!isinteger(BASE)) { push(BASE); numerator(); push(EXPO); power_numbers(); push(BASE); denominator(); push(EXPO); negate(); power_numbers(); multiply(); return; } // At this point BASE is a positive integer. // If BASE is small then factor it. if (is_small_integer(BASE)) { push(BASE); push(EXPO); quickfactor(); return; } // At this point BASE is a positive integer and EXPO is not an integer. if (MLENGTH(EXPO->u.q.a) > 1 || MLENGTH(EXPO->u.q.b) > 1) { push_symbol(POWER); push(BASE); push(EXPO); list(3); return; } a = EXPO->u.q.a[0]; b = EXPO->u.q.b[0]; x = mroot(BASE->u.q.a, b); if (x == 0) { push_symbol(POWER); push(BASE); push(EXPO); list(3); return; } y = mpow(x, a); mfree(x); p3 = alloc(); p3->k = NUM; if (MSIGN(EXPO->u.q.a) == -1) { p3->u.q.a = mint(1); p3->u.q.b = y; } else { p3->u.q.a = y; p3->u.q.b = mint(1); } push(p3); } //----------------------------------------------------------------------------- // // Normalize the angle of unit imaginary, i.e. (-1) ^ N // // Input: N on stack (must be rational, not float) // // Output: Result on stack // // Note: // // n = q * d + r // // Example: // n d q r // // (-1)^(8/3) -> (-1)^(2/3) 8 3 2 2 // (-1)^(7/3) -> (-1)^(1/3) 7 3 2 1 // (-1)^(5/3) -> -(-1)^(2/3) 5 3 1 2 // (-1)^(4/3) -> -(-1)^(1/3) 4 3 1 1 // (-1)^(2/3) -> (-1)^(2/3) 2 3 0 2 // (-1)^(1/3) -> (-1)^(1/3) 1 3 0 1 // // (-1)^(-1/3) -> -(-1)^(2/3) -1 3 -1 2 // (-1)^(-2/3) -> -(-1)^(1/3) -2 3 -1 1 // (-1)^(-4/3) -> (-1)^(2/3) -4 3 -2 2 // (-1)^(-5/3) -> (-1)^(1/3) -5 3 -2 1 // (-1)^(-7/3) -> -(-1)^(2/3) -7 3 -3 2 // (-1)^(-8/3) -> -(-1)^(1/3) -8 3 -3 1 // //----------------------------------------------------------------------------- #define A p1 #define Q p2 #define R p3 static void normalize_angle(void) { save(); A = pop(); // integer exponent? if (isinteger(A)) { if (A->u.q.a[0] & 1) push_integer(-1); // odd exponent else push_integer(1); // even exponent restore(); return; } // floor push(A); bignum_truncate(); Q = pop(); if (isnegativenumber(A)) { push(Q); push_integer(-1); add(); Q = pop(); } // remainder (always positive) push(A); push(Q); subtract(); R = pop(); // remainder becomes new angle push_symbol(POWER); push_integer(-1); push(R); list(3); // negate if quotient is odd if (Q->u.q.a[0] & 1) negate(); restore(); } static int is_small_integer(U *p) { if (isinteger(p) && MLENGTH(p->u.q.a) == 1 && (p->u.q.a[0] & 0x80000000) == 0) return 1; else return 0; }