Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
Implement roots() builtins
Signed-off-by: Moritz Lenz <moritz@faui2k3.org>
  • Loading branch information
leto authored and moritz committed Mar 24, 2009
1 parent 07ceac6 commit 1ef8919
Show file tree
Hide file tree
Showing 3 changed files with 116 additions and 0 deletions.
7 changes: 7 additions & 0 deletions CREDITS
Expand Up @@ -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
Expand Down
90 changes: 90 additions & 0 deletions src/builtins/any-num.pir
Expand Up @@ -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
Expand Down
19 changes: 19 additions & 0 deletions src/builtins/math.pir
Expand Up @@ -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 )
Expand All @@ -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


Expand All @@ -62,6 +79,8 @@ constant I<e>.

&log10 := &log.assuming:base(10);

Returns the base 10 logarithm of $x.

=cut

.sub 'log10'
Expand Down

0 comments on commit 1ef8919

Please sign in to comment.