nisaguci.pkl
/
aljabar
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
aljabar/node_modules/algebrite/sources/multiply.coffee

620 lines
13 KiB
CoffeeScript
Raw Permalink Blame History

# Symbolic multiplication
# multiplication is commutative, so it can't be used
# e.g. on two matrices.
# But it can be used, say, on a scalar and a matrix.,
# so the output of a multiplication is not
# always a scalar.
#extern void append(void)
#static void parse_p1(void)
#static void parse_p2(void)
#static void __normalize_radical_factors(int)
DEBUG_MULTIPLY = false
Eval_multiply = ->
push(cadr(p1))
Eval()
p1 = cddr(p1)
while (iscons(p1))
push(car(p1))
Eval()
multiply()
p1 = cdr(p1)
# this one doesn't eval the factors,
# so you pass i*(-1)^(1/2), it wouldnt't
# give -1, because i is not evalled
multiply = ->
if (esc_flag)
stop("escape key stop")
if (isNumericAtom(stack[tos - 2]) && isNumericAtom(stack[tos - 1]))
multiply_numbers()
else
save()
yymultiply()
restore()
yymultiply = ->
h = 0
i = 0
n = 0
# pop operands
p2 = pop()
p1 = pop()
h = tos
# is either operand zero?
if (isZeroAtom(p1) || isZeroAtom(p2))
if evaluatingAsFloats then push_double(0.0) else push(zero)
return
# is either operand a sum?
#console.log("yymultiply: expanding: " + expanding)
if (expanding && isadd(p1))
p1 = cdr(p1)
if evaluatingAsFloats then push_double(0.0) else push(zero)
while (iscons(p1))
push(car(p1))
push(p2)
multiply()
add()
p1 = cdr(p1)
return
if (expanding && isadd(p2))
p2 = cdr(p2)
if evaluatingAsFloats then push_double(0.0) else push(zero)
while (iscons(p2))
push(p1)
push(car(p2))
multiply()
add()
p2 = cdr(p2)
return
# scalar times tensor?
if (!istensor(p1) && istensor(p2))
push(p1)
push(p2)
scalar_times_tensor()
return
# tensor times scalar?
if (istensor(p1) && !istensor(p2))
push(p1)
push(p2)
tensor_times_scalar()
return
# adjust operands
if (car(p1) == symbol(MULTIPLY))
p1 = cdr(p1)
else
push(p1)
list(1)
p1 = pop()
if (car(p2) == symbol(MULTIPLY))
p2 = cdr(p2)
else
push(p2)
list(1)
p2 = pop()
# handle numerical coefficients
if (isNumericAtom(car(p1)) && isNumericAtom(car(p2)))
push(car(p1))
push(car(p2))
multiply_numbers()
p1 = cdr(p1)
p2 = cdr(p2)
else if (isNumericAtom(car(p1)))
push(car(p1))
p1 = cdr(p1)
else if (isNumericAtom(car(p2)))
push(car(p2))
p2 = cdr(p2)
else
if evaluatingAsFloats then push_double(1.0) else push(one)
parse_p1()
parse_p2()
while (iscons(p1) && iscons(p2))
# if (car(p1)->gamma && car(p2)->gamma) {
# combine_gammas(h)
# p1 = cdr(p1)
# p2 = cdr(p2)
# parse_p1()
# parse_p2()
# continue
# }
if (caar(p1) == symbol(OPERATOR) && caar(p2) == symbol(OPERATOR))
push_symbol(OPERATOR)
push(cdar(p1))
push(cdar(p2))
append()
cons()
p1 = cdr(p1)
p2 = cdr(p2)
parse_p1()
parse_p2()
continue
switch (cmp_expr(p3, p4))
when -1
push(car(p1))
p1 = cdr(p1)
parse_p1()
when 1
push(car(p2))
p2 = cdr(p2)
parse_p2()
when 0
combine_factors(h)
p1 = cdr(p1)
p2 = cdr(p2)
parse_p1()
parse_p2()
else
stop("internal error 2")
# push remaining factors, if any
while (iscons(p1))
push(car(p1))
p1 = cdr(p1)
while (iscons(p2))
push(car(p2))
p2 = cdr(p2)
# normalize radical factors
# example: 2*2(-1/2) -> 2^(1/2)
# must be done after merge because merge may produce radical
# example: 2^(1/2-a)*2^a -> 2^(1/2)
__normalize_radical_factors(h)
# this hack should not be necessary, unless power returns a multiply
#for (i = h; i < tos; i++) {
# if (car(stack[i]) == symbol(MULTIPLY)) {
# multiply_all(tos - h)
# return
# }
#}
if (expanding)
for i in [h...tos]
if (isadd(stack[i]))
multiply_all(tos - h)
return
# n is the number of result factors on the stack
n = tos - h
if (n == 1)
return
# discard integer 1
if (isrational(stack[h]) && equaln(stack[h], 1))
if (n == 2)
p7 = pop()
pop()
push(p7)
else
stack[h] = symbol(MULTIPLY)
list(n)
return
list(n)
p7 = pop()
push_symbol(MULTIPLY)
push(p7)
cons()
# Decompose a factor into base and power.
#
# input: car(p1) factor
#
# output: p3 factor's base
#
# p5 factor's power (possibly 1)
parse_p1 = ->
p3 = car(p1)
p5 = if evaluatingAsFloats then one_as_double else one
if (car(p3) == symbol(POWER))
p5 = caddr(p3)
p3 = cadr(p3)
# Decompose a factor into base and power.
#
# input: car(p2) factor
#
# output: p4 factor's base
#
# p6 factor's power (possibly 1)
parse_p2 = ->
p4 = car(p2)
p6 = if evaluatingAsFloats then one_as_double else one
if (car(p4) == symbol(POWER))
p6 = caddr(p4)
p4 = cadr(p4)
# h an integer
combine_factors = (h) ->
push(p4)
push(p5)
push(p6)
add()
power()
p7 = pop()
if (isNumericAtom(p7))
push(stack[h])
push(p7)
multiply_numbers()
stack[h] = pop()
else if (car(p7) == symbol(MULTIPLY))
# power can return number * factor (i.e. -1 * i)
if (isNumericAtom(cadr(p7)) && cdddr(p7) == symbol(NIL))
push(stack[h])
push(cadr(p7))
multiply_numbers()
stack[h] = pop()
push(caddr(p7))
else
push(p7)
else
push(p7)
gp = [
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,1,-6,-7,-8,-3,-4,-5,13,14,15,-16,9,10,11,-12],
[0,0,6,-1,-11,10,-2,-15,14,12,-5,4,-9,16,-8,7,-13],
[0,0,7,11,-1,-9,15,-2,-13,5,12,-3,-10,8,16,-6,-14],
[0,0,8,-10,9,-1,-14,13,-2,-4,3,12,-11,-7,6,16,-15],
[0,0,3,2,15,-14,1,11,-10,16,-8,7,13,12,-5,4,9],
[0,0,4,-15,2,13,-11,1,9,8,16,-6,14,5,12,-3,10],
[0,0,5,14,-13,2,10,-9,1,-7,6,16,15,-4,3,12,11],
[0,0,13,12,-5,4,16,-8,7,-1,-11,10,-3,-2,-15,14,-6],
[0,0,14,5,12,-3,8,16,-6,11,-1,-9,-4,15,-2,-13,-7],
[0,0,15,-4,3,12,-7,6,16,-10,9,-1,-5,-14,13,-2,-8],
[0,0,16,-9,-10,-11,-13,-14,-15,-3,-4,-5,1,-6,-7,-8,2],
[0,0,9,-16,8,-7,-12,5,-4,-2,-15,14,6,-1,-11,10,3],
[0,0,10,-8,-16,6,-5,-12,3,15,-2,-13,7,11,-1,-9,4],
[0,0,11,7,-6,-16,4,-3,-12,-14,13,-2,8,-10,9,-1,5],
[0,0,12,13,14,15,9,10,11,-6,-7,-8,-2,-3,-4,-5,-1]
]
#if 0
# h an int
combine_gammas = (h) ->
n = gp[Math.floor(p1.gamma)][Math.floor(p2.gamma)]
if (n < 0)
n = -n
push(stack[h])
negate()
stack[h] = pop()
if (n > 1)
push(_gamma[n])
# this is useful for example when you are just adding/removing
# factors from an already factored quantity.
# e.g. if you factored x^2 + 3x + 2 into (x+1)(x+2)
# and you want to divide by (x+1) , i.e. you multiply by (x-1)^-1,
# then there is no need to expand.
multiply_noexpand = ->
prev_expanding = expanding
expanding = 0
multiply()
expanding = prev_expanding
# multiply n factors on stack
# n an integer
multiply_all = (n) ->
i = 0
if (n == 1)
return
if (n == 0)
push if evaluatingAsFloats then one_as_double else one
return
h = tos - n
push(stack[h])
for i in [1...n]
push(stack[h + i])
multiply()
stack[h] = pop()
moveTos h + 1
# n an integer
multiply_all_noexpand = (n) ->
prev_expanding = expanding
expanding = 0
multiply_all(n)
expanding = prev_expanding
#-----------------------------------------------------------------------------
#
# Symbolic division, or numeric division if doubles are found.
#
# Input: Dividend and divisor on stack
#
# Output: Quotient on stack
#
#-----------------------------------------------------------------------------
divide = ->
if (isNumericAtom(stack[tos - 2]) && isNumericAtom(stack[tos - 1]))
divide_numbers()
else
inverse()
multiply()
# this is different from inverse of a matrix (inv)!
inverse = ->
if (isNumericAtom(stack[tos - 1]))
invert_number()
else
push_integer(-1)
power()
reciprocate = ->
inverse()
negate = ->
if (isNumericAtom(stack[tos - 1]))
negate_number()
else
if evaluatingAsFloats then push_double(-1.0) else push_integer(-1)
multiply()
negate_expand = ->
prev_expanding = expanding
expanding = 1
negate()
expanding = prev_expanding
negate_noexpand = ->
prev_expanding = expanding
expanding = 0
negate()
expanding = prev_expanding
#-----------------------------------------------------------------------------
#
# Normalize radical factors
#
# Input: stack[h] Coefficient factor, possibly 1
#
# stack[h + 1] Second factor
#
# stack[tos - 1] Last factor
#
# Output: Reduced coefficent and normalized radicals (maybe)
#
# Example: 2*2^(-1/2) -> 2^(1/2)
#
# (power number number) is guaranteed to have the following properties:
#
# 1. Base is an integer
#
# 2. Absolute value of exponent < 1
#
# These properties are assured by the power function.
#
#-----------------------------------------------------------------------------
#define A p1
#define B p2
#define BASE p3
#define EXPO p4
#define TMP p5
# h is an int
__normalize_radical_factors = (h) ->
i = 0
# if coeff is 1 or floating then don't bother
if (isplusone(stack[h]) || isminusone(stack[h]) || isdouble(stack[h]))
return
# if no radicals then don't bother
for i in [(h + 1)...tos]
if (__is_radical_number(stack[i]))
break
if (i == tos)
return
# ok, try to simplify
save()
# numerator
push(stack[h])
mp_numerator()
if DEBUG_MULTIPLY then console.log("__normalize_radical_factors numerator: " + stack[tos-1])
p1 = pop(); # p1 is A
for i in [(h + 1)...tos]
if (isplusone(p1) || isminusone(p1)) # p1 is A
break
if (!__is_radical_number(stack[i]))
continue
p3 = cadr(stack[i]); #p3 is BASE
p4 = caddr(stack[i]); #p4 is EXPO
# exponent must be negative
if (!isnegativenumber(p4)) #p4 is EXPO
continue
# numerator divisible by p3 (base)?
push(p1); # p1 is A
push(p3); #p3 is BASE
divide()
p5 = pop(); #p5 is TMP
if (!isinteger(p5)) #p5 is TMP
continue
# reduce numerator
p1 = p5; # p1 is A #p5 is TMP
# invert radical
push_symbol(POWER)
push(p3); #p3 is BASE
push if evaluatingAsFloats then one_as_double else one
push(p4); #p4 is EXPO
add()
list(3)
stack[i] = pop()
# denominator
push(stack[h])
mp_denominator()
if DEBUG_MULTIPLY then console.log("__normalize_radical_factors denominator: " + stack[tos-1])
p2 = pop(); # p2 is B
for i in [(h + 1)...tos]
if (isplusone(p2)) # p2 is B
break
if (!__is_radical_number(stack[i]))
continue
p3 = cadr(stack[i]); #p3 is BASE
p4 = caddr(stack[i]); #p4 is EXPO
# exponent must be positive
if (isnegativenumber(p4)) #p4 is EXPO
continue
# denominator divisible by p3? #p3 is BASE
push(p2); # p2 is B
push(p3); #p3 is BASE
divide()
p5 = pop(); #p5 is TMP
if (!isinteger(p5)) #p5 is TMP
continue
if DEBUG_MULTIPLY then console.log("__new radical p5: " + p5.toString())
if DEBUG_MULTIPLY then console.log("__new radical top stack: " + stack[tos-1])
# reduce denominator
p2 = p5; # p2 is B #p5 is TMP
# invert radical
push_symbol(POWER)
push(p3); #p3 is BASE
push(p4); #p4 is EXPO
if DEBUG_MULTIPLY then console.log("__new radical p3: " + p3.toString())
if DEBUG_MULTIPLY then console.log("__new radical p4: " + p4.toString())
push(one)
subtract()
if dontCreateNewRadicalsInDenominatorWhenEvalingMultiplication
if (isinteger(p3) and !isinteger(stack[tos-1]) and isnegativenumber(stack[tos - 1]))
# bail out,
# we want to avoid going ahead with the subtraction of
# the exponents, because that would turn a perfectly good
# integer exponent in the denominator into a fractional one
# i.e. a radical.
# Note that this only prevents new radicals ending up
# in the denominator, it doesn't fix existing ones.
pop()
pop()
pop()
push(p1); # p1 is A
push(p3); #p3 is BASE
divide()
p1 = pop()
break
if DEBUG_MULTIPLY then console.log("__new radical exponent: " + stack[tos-1])
list(3)
stack[i] = pop()
# reconstitute the coefficient
push(p1); # p1 is A
push(p2); # p2 is B
divide()
stack[h] = pop()
restore()
# don't include i
# p is a U
# TODO should this be in is.coffee ?
__is_radical_number = (p) ->
# don't use i
car(p) == symbol(POWER) && isNumericAtom(cadr(p)) && isfraction(caddr(p)) && !isminusone(cadr(p))
#-----------------------------------------------------------------------------
#
# > a*hilbert(2)
# ((a,1/2*a),(1/2*a,1/3*a))
#
# Note that "a" is presumed to be a scalar. Is this correct?
#
# Yes, because "*" has no meaning if "a" is a tensor.
# To multiply tensors, "dot" or "outer" should be used.
#
# > dot(a,hilbert(2))
# dot(a,((1,1/2),(1/2,1/3)))
#
# In this case "a" could be a scalar or tensor so the result is not
# expanded.
#
#-----------------------------------------------------------------------------
Reference in New Issue View Git Blame Copy Permalink