Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

Implement roots() builtins

Signed-off-by: Moritz Lenz <moritz@faui2k3.org>
  • Loading branch information...
commit 1ef891980c196cd217f12de6b02a4c1e963993d1 1 parent 07ceac6
@leto leto authored moritz committed
Showing with 116 additions and 0 deletions.
  1. +7 −0 CREDITS
  2. +90 −0 src/builtins/any-num.pir
  3. +19 −0 src/builtins/math.pir
View
7 CREDITS
@@ -181,6 +181,13 @@ N: Jonathan Scott Duff
U: duff
E: duff@pobox.com
+N: Jonathan "Duke" Leto
+U: leto
+U: dukeleto
+E: jonathan@leto.net
+W: http://leto.net
+S: Portland, OR
+
N: Jonathan Worthington
U: jonathan
E: jnthn@jnthn.net
View
90 src/builtins/any-num.pir
@@ -184,6 +184,96 @@ error.
=cut
+=item roots
+
+ our Array multi Num::roots ( Complex $z, Int $n )
+
+Returns an Array consisting of the $n roots of a Complex number $z, where $n is
+a positive integer. For any Complex number $z ( which includes real numbers and
+integers as a subset ) there are a set of $n numbers W such that $w_k ** $n = $z,
+or in set theory notation:
+
+W = { $w_i : $w_i ** $n = $z and 0 <= i <= n-1 } .
+
+These can be written in terms of the multiple-valued complex logarithm:
+
+exp[1/$n*(Log($z)]
+
+which is equal to
+
+$w_k = exp[1/$n*(log($r)+i*($theta + 2*k*pi))] where k = 0,1,2,..., n-1
+
+where ($r,$theta) = $z.polar . The angle $theta returned is always in the
+interval -pi <= $theta <= pi .
+
+=cut
+
+.sub 'roots' :method
+ .param int n
+ .local num pi, r, theta
+ .local pmc x, result, roots
+ x = self
+ pi = atan 1
+ pi *= 4
+ roots = new 'FixedPMCArray'
+ if n > 0 goto positive
+ roots = 1 # single element array
+ roots[0] = 'NaN'
+ goto done
+ positive:
+ roots = n # fix array size to n
+ if n > 1 goto general
+ roots[0] = x
+ goto done
+ general:
+ div $N0, 1, n
+ $I0 = 0
+ $I1 = isa x, 'Complex'
+ unless $I1 goto real
+ $N6 = x[0]
+ $N7 = x[1]
+ theta = atan $N7, $N6 # angle of polar representation
+ $N6 *= $N6
+ $N7 *= $N7
+ $N8 = $N6 + $N7
+ r = sqrt $N8 # radius of polar representation
+ $N1 = ln r
+ goto loop
+ real:
+ $N4 = x
+ $N4 = abs $N4 # if x < 0 we rotate by exp(i pi/n) later on
+ $N1 = ln $N4 # ln(abs(x)) = ln(r)
+ loop:
+ if $I0 >= n goto done
+ $P2 = new 'Complex' # this can surely be optimized
+ $N3 = $N0
+ $N3 *= 2
+ $N3 *= pi
+ $N3 *= $I0
+ $P2[1] = $N3 # 2*$I0*pi/n
+ $N5 = $P2[1]
+ unless $I1 goto rotate_negative_reals
+ $N8 = $N0
+ $N8 *= theta # theta/n
+ $N5 += $N8 # 2*$I0*pi/n + theta/n
+ goto exponentiate
+ rotate_negative_reals: # we must rotate answer since we factored out (-1)^(1/n)
+ if x > 0 goto exponentiate
+ div $N4, pi, n
+ $N5 += $N4 # exp( i pi / n ) = (-1)^(1/n) (principle root)
+ exponentiate:
+ $N9 = $N0
+ $N9 *= $N1 # 1/n*ln(r)
+ $P2[0] = $N9
+ $P2[1] = $N5
+ $P2 = $P2.'exp'() # exp(1/n*(ln(r)+i*(theta + 2*k*pi)))
+ roots[$I0] = $P2
+ inc $I0
+ goto loop
+ done:
+ .return (roots)
+.end
+
.sub 'unpolar' :method
.param num angle
.local num mag
View
19 src/builtins/math.pir
@@ -15,6 +15,15 @@ src/builtins/math.pir - Perl6 math functions
## TODO: figure out what to get working, in order to uncomment the following
## .namespace [ 'Math::Basic' ]
+
+.sub 'roots' :multi(_, 'Integer')
+ .param pmc x
+ .param int n
+ .local pmc result
+ result = x.'roots'(n)
+ .return (result)
+.end
+
=item sign
our Int multi Num::sign ( Num $x )
@@ -32,12 +41,20 @@ or more succinctly:
$x <=> 0;
}
+Returns the sign of $x, i.e +1 for positive numbers (including Inf), zero for zero and -1 for negative numbers (including -Inf).
+
=cut
.sub 'sign'
.param pmc a
+ if a == 'Inf' goto unity
+ if a == 'NaN' goto not_a_number
$I0 = cmp_num a, 0
.return ($I0)
+ not_a_number:
+ .return (a)
+ unity:
+ .return (1)
.end
@@ -62,6 +79,8 @@ constant I<e>.
&log10 := &log.assuming:base(10);
+Returns the base 10 logarithm of $x.
+
=cut
.sub 'log10'
Please sign in to comment.
Something went wrong with that request. Please try again.