Roots of unity
You are encouraged to solve this task according to the task description, using any language you may know.
The purpose of this task is to explore working with complex numbers. Given n, find the n-th roots of unity.
Ada
<lang ada> with Ada.Text_IO; use Ada.Text_IO; with Ada.Float_Text_IO; use Ada.Float_Text_IO; with Ada.Numerics.Complex_Types; use Ada.Numerics.Complex_Types;
procedure Roots_Of_Unity is
Root : Complex;
begin
for N in 2..10 loop Put_Line ("N =" & Integer'Image (N)); for K in 0..N - 1 loop Root := Compose_From_Polar ( Modulus => 1.0, Argument => Float (K), Cycle => Float (N) ); -- Output Put (" k =" & Integer'Image (K) & ", "); if Re (Root) < 0.0 then Put ("-"); else Put ("+"); end if; Put (abs Re (Root), Fore => 1, Exp => 0); if Im (Root) < 0.0 then Put ("-"); else Put ("+"); end if; Put (abs Im (Root), Fore => 1, Exp => 0); Put_Line ("i"); end loop; end loop;
end Roots_Of_Unity; </lang> Ada provides a direct implementation of polar composition of complex numbers x e2πi y. The function Compose_From_Polar is used to compose roots. The third argument of the function is the cycle. Instead of the standard cycle 2π, N is used. Sample output:
N = 2 k = 0, +1.00000+0.00000i k = 1, -1.00000+0.00000i N = 3 k = 0, +1.00000+0.00000i k = 1, -0.50000+0.86603i k = 2, -0.50000-0.86603i N = 4 k = 0, +1.00000+0.00000i k = 1, +0.00000+1.00000i k = 2, -1.00000+0.00000i k = 3, +0.00000-1.00000i N = 5 k = 0, +1.00000+0.00000i k = 1, +0.30902+0.95106i k = 2, -0.80902+0.58779i k = 3, -0.80902-0.58779i k = 4, +0.30902-0.95106i N = 6 k = 0, +1.00000+0.00000i k = 1, +0.50000+0.86603i k = 2, -0.50000+0.86603i k = 3, -1.00000+0.00000i k = 4, -0.50000-0.86603i k = 5, +0.50000-0.86603i N = 7 k = 0, +1.00000+0.00000i k = 1, +0.62349+0.78183i k = 2, -0.22252+0.97493i k = 3, -0.90097+0.43388i k = 4, -0.90097-0.43388i k = 5, -0.22252-0.97493i k = 6, +0.62349-0.78183i N = 8 k = 0, +1.00000+0.00000i k = 1, +0.70711+0.70711i k = 2, +0.00000+1.00000i k = 3, -0.70711+0.70711i k = 4, -1.00000+0.00000i k = 5, -0.70711-0.70711i k = 6, +0.00000-1.00000i k = 7, +0.70711-0.70711i N = 9 k = 0, +1.00000+0.00000i k = 1, +0.76604+0.64279i k = 2, +0.17365+0.98481i k = 3, -0.50000+0.86603i k = 4, -0.93969+0.34202i k = 5, -0.93969-0.34202i k = 6, -0.50000-0.86603i k = 7, +0.17365-0.98481i k = 8, +0.76604-0.64279i N = 10 k = 0, +1.00000+0.00000i k = 1, +0.80902+0.58779i k = 2, +0.30902+0.95106i k = 3, -0.30902+0.95106i k = 4, -0.80902+0.58779i k = 5, -1.00000+0.00000i k = 6, -0.80902-0.58779i k = 7, -0.30902-0.95106i k = 8, +0.30902-0.95106i k = 9, +0.80902-0.58779i
ALGOL 68
FORMAT complex fmt=$g(-6,4)"⊥"g(-6,4)$; FOR root FROM 2 TO 10 DO printf(($g(4)$,root)); FOR n TO root OVER 2 DO printf(($xf(complex fmt)$,complex exp( 0 I 2*pi*n/root))) OD; printf($l$) OD
Output:
+2 -1.000⊥0.0000 +3 -.5000⊥0.8660 +4 0.0000⊥1.0000 -1.000⊥0.0000 +5 0.3090⊥0.9511 -.8090⊥0.5878 +6 0.5000⊥0.8660 -.5000⊥0.8660 -1.000⊥0.0000 +7 0.6235⊥0.7818 -.2225⊥0.9749 -.9010⊥0.4339 +8 0.7071⊥0.7071 0.0000⊥1.0000 -.7071⊥0.7071 -1.000⊥0.0000 +9 0.7660⊥0.6428 0.1736⊥0.9848 -.5000⊥0.8660 -.9397⊥0.3420 +10 0.8090⊥0.5878 0.3090⊥0.9511 -.3090⊥0.9511 -.8090⊥0.5878 -1.000⊥0.0000
BASIC
For high n's, this may repeat the root of 1 + 0*i.
CLS PI = 3.1415926# n = 5 'this can be changed for any desired n angle = 0 'start at angle 0 DO real = COS(angle) 'real axis is the x axis IF (ABS(real) < 10 ^ -5) THEN real = 0 'get rid of annoying sci notation imag = SIN(angle) 'imaginary axis is the y axis IF (ABS(imag) < 10 ^ -5) THEN imag = 0 'get rid of annoying sci notation PRINT real; "+"; imag; "i" 'answer on every line angle = angle + (2 * PI) / n 'all the way around the circle at even intervals LOOP WHILE angle < 2 * PI
C
<lang c> #include <stdio.h>
#include <math.h> #define PI 3.1415926 int main (int argc, char *argv[]) { char sign; int i, n; float rpart, ipart, angle; for (n = 2; n <= 10; n++) { angle = 0.0; printf("%d: ", n); for (i = 1; i <= n; i++) { rpart = cos(angle); ipart = sin(angle); if (ipart < 0) sign = '-'; else sign = '+'; printf("%5.4f%cj%5.4f ", rpart, sign, fabs(ipart)); angle = angle + 2.0*PI/(float)n; } printf("\n"); } }
</lang>
C++
<lang cpp>#include <complex>
- include <cmath>
- include <iostream>
double const pi = 4 * std::atan(1);
int main() {
for (int n = 2; n <= 10; ++n) { std::cout << n << ": "; for (int k = 0; k < n; ++k) std::cout << std::polar(1, 2*pi*k/n) << " "; std::cout << std::endl; }
}</lang>
Common Lisp
<lang lisp>(defun roots-of-unity (n)
(loop for i below n collect (cis (* pi (/ (* 2 i) n)))))</lang>
The expression is slightly more complicated than necessary in order to preserve exact rational arithmetic until multiplying by pi. The author of this example is not a floating point expert and not sure whether this is actually useful; if not, the simpler expression is (cis (/ (* 2 pi i) n)).
D
<lang d>module nthroots ; import std.stdio, std.math ;
creal[] nthroots(int n) {
creal[] res ; for(int k = 1 ; k <= n ; k++) res ~= expi(PI*2*k/n) ; return res ;
} void main() {
for(int i = 1; i <= 8 ; i++) writefln("%2dth : %5.2f", i, nthroots(i)) ;
}</lang>
Forth
Complex numbers are not a native type in Forth, so we calculate the roots by hand.
: f0. ( f -- ) fdup 0e 0.001e f~ if fdrop 0e then f. ; : .roots ( n -- ) dup 1 do pi i 2* 0 d>f f* dup 0 d>f f/ ( F: radians ) fsincos cr ." real " f0. ." imag " f0. loop drop ; 3 set-precision 5 .roots
Fortran
Sin/Cos + Scalar Loop
<lang fortran>
PROGRAM Roots COMPLEX :: root INTEGER :: i, n REAL :: angle, pi pi = 4.0 * ATAN(1.0) DO n = 2, 7 angle = 0.0 WRITE(*,"(I1,A)", ADVANCE="NO") n,": " DO i = 1, n root = CMPLX(COS(angle), SIN(angle)) WRITE(*,"(SP,2F7.4,A)", ADVANCE="NO") root, "j " angle = angle + (2.0*pi / REAL(n)) END DO WRITE(*,*) END DO END PROGRAM Roots
</lang> Output
2: +1.0000+0.0000j -1.0000+0.0000j 3: +1.0000+0.0000j -0.5000+0.8660j -0.5000-0.8660j 4: +1.0000+0.0000j +0.0000+1.0000j -1.0000+0.0000j +0.0000-1.0000j 5: +1.0000+0.0000j +0.3090+0.9511j -0.8090+0.5878j -0.8090-0.5878j +0.3090-0.9511j 6: +1.0000+0.0000j +0.5000+0.8660j -0.5000+0.8660j -1.0000+0.0000j -0.5000-0.8660j +0.5000-0.8660j 7: +1.0000+0.0000j +0.6235+0.7818j -0.2225+0.9749j -0.9010+0.4339j -0.9010-0.4339j -0.2225-0.9749j +0.6235-0.7818j
Exp + Array-valued Statement
<lang fortran> program unity
real, parameter :: pi = 3.141592653589793 complex, parameter :: i = (0, 1) complex, dimension(0:7-1) :: unit_circle integer :: n, j do n = 2, 7 !!!! KEY STEP, does all the calculations in one statement !!!! unit_circle(0:n-1) = exp(2*i*pi/n * (/ (j, j=0, n-1) /) )
write(*,"(i1,a)", advance="no") n, ": " write(*,"(sp,2f7.4,a)", advance="no") (unit_circle(j), "j ", j = 0, n-1) write(*,*) end do end program unity</lang>
Haskell
import Data.Complex rootsOfUnity n = [mkPolar 1.0 (2*pi*k/n) | k <- [1..n]]
Output:
*Main> rootsOfUnity 3 [(-0.4999999999999998) :+ 0.8660254037844387, (-0.5000000000000004) :+ (-0.8660254037844384), 1.0 :+ (-2.4492127076447545e-16)]
Icon
While icon does have complex type, it does not know about polar
procedure main() roots(10) end procedure roots(n) every n := 2 to 10 do every writes(n | (str_rep((0 to (n-1)) * 2 * &pi / n)) | "\n") end procedure str_rep(k) return " " || cos(k) || "+" || sin(k) || "i" end
IDL
For some example n:
n = 5 print, exp( dcomplex( 0, 2*!pi/n) ) ^ ( 1 + indgen(n) )
Outputs:
( 0.30901696, 0.95105653)( -0.80901704, 0.58778520)( -0.80901693, -0.58778534)( 0.30901713, -0.95105647)( 1.0000000, 1.7484556e-007)
J
rou=: [: ^ i. * (o. 0j2) % ] rou 4 1 0j1 _1 0j_1 rou 5 1 0.309017j0.951057 _0.809017j0.587785 _0.809017j_0.587785 0.309017j_0.951057
The computation can also be written as a loop, shown here for comparison only.
rou1=: 3 : 0 z=. 0 $ r=. ^ o. 0j2 % y [ e=. 1 for. i.y do. z=. z,e e=. e*r end. z )
Java
Java doesn't have a nice way of dealing with complex numbers, so the real and imaginary parts are calculated separately based on the angle and printed together. There are also checks in this implementation to get rid of extremely small values (< 1.0E-3 where scientific notation sets in for Doubles). Instead, they are simply represented as 0. To remove those checks (for very high n's), remove both if statements. <lang java>public static void unity(int n){ //all the way around the circle at even intervals for(double angle = 0;angle < 2 * Math.PI;angle += (2 * Math.PI) / n){ double real = Math.cos(angle); //real axis is the x axis if(Math.abs(real) < 1.0E-3) real = 0.0; //get rid of annoying sci notation double imag = Math.sin(angle); //imaginary axis is the y axis if(Math.abs(imag) < 1.0E-3) imag = 0.0; //get rid of annoying sci notation System.out.print(real + " + " + imag + "i\t"); //tab-separated answers } }</lang>
OCaml
<lang ocaml>open Complex
let pi = 4. *. atan 1.
let () =
for n = 1 to 10 do Printf.printf "%2d " n; for k = 1 to n do let ret = polar 1. (2. *. pi *. float_of_int k /. float_of_int n) in Printf.printf "(%f + %f i)" ret.re ret.im done; print_newline () done</lang>
Octave
<lang octave>for j = 2 : 10
printf("*** %d\n", j); for n = 1 : j disp(exp(2i*pi*n/j)); endfor disp("");
endfor</lang>
Perl
<lang perl>use Math::Complex;
foreach $n (2 .. 10) {
printf "%2d", $n; foreach $k (0 .. $n-1) { $ret = cplxe(1, 2 * pi * $k / $n); $ret->display_format('style' => 'cartesian', 'format' => '%.3f'); print " $ret"; } print "\n";
}</lang> Output:
2 1.000 -1.000+0.000i 3 1.000 -0.500+0.866i -0.500-0.866i 4 1.000 0.000+1.000i -1.000+0.000i -0.000-1.000i 5 1.000 0.309+0.951i -0.809+0.588i -0.809-0.588i 0.309-0.951i 6 1.000 0.500+0.866i -0.500+0.866i -1.000+0.000i -0.500-0.866i 0.500-0.866i 7 1.000 0.623+0.782i -0.223+0.975i -0.901+0.434i -0.901-0.434i -0.223-0.975i 0.623-0.782i 8 1.000 0.707+0.707i 0.000+1.000i -0.707+0.707i -1.000+0.000i -0.707-0.707i -0.000-1.000i 0.707-0.707i 9 1.000 0.766+0.643i 0.174+0.985i -0.500+0.866i -0.940+0.342i -0.940-0.342i -0.500-0.866i 0.174-0.985i 0.766-0.643i 10 1.000 0.809+0.588i 0.309+0.951i -0.309+0.951i -0.809+0.588i -1.000+0.000i -0.809-0.588i -0.309-0.951i 0.309-0.951i 0.809-0.588i
Python
<lang python> import cmath class Complex(complex):
def __repr__(self): rp = '%7.5f'%self.real if not self.pureImag() else ip = '%7.5fj'%self.imag if not self.pureReal() else conj = if (self.pureImag() or self.pureReal() or self.imag<0.0) else '+' return '0.0' if (self.pureImag() and self.pureReal()) else rp+conj+ip def pureImag(self): return abs( self.real) < 0.000005 def pureReal(self): return abs( self.imag) < 0.000005
def croots(n):
if n<=0: return None return (Complex(cmath.exp(2j*k*cmath.pi/n)) for k in range(n)) # in Python 2.6+: return (Complex(cmath.rect(1, 2*k*cmath.pi/n)) for k in range(n))
for nr in range(2,11):
print nr, list(croots(nr))
</lang> Output:
2 [1.00000, -1.00000] 3 [1.00000, -0.50000+0.86603j, -0.50000-0.86603j] 4 [1.00000, 1.00000j, -1.00000, -1.00000j] 5 [1.00000, 0.30902+0.95106j, -0.80902+0.58779j, -0.80902-0.58779j, 0.30902-0.95106j] 6 [1.00000, 0.50000+0.86603j, -0.50000+0.86603j, -1.00000, -0.50000-0.86603j, 0.50000-0.86603j] 7 [1.00000, 0.62349+0.78183j, -0.22252+0.97493j, -0.90097+0.43388j, -0.90097-0.43388j, -0.22252-0.97493j, 0.62349-0.78183j] 8 [1.00000, 0.70711+0.70711j, 1.00000j, -0.70711+0.70711j, -1.00000, -0.70711-0.70711j, -1.00000j, 0.70711-0.70711j] 9 [1.00000, 0.76604+0.64279j, 0.17365+0.98481j, -0.50000+0.86603j, -0.93969+0.34202j, -0.93969-0.34202j, -0.50000-0.86603j, 0.17365-0.98481j, 0.76604-0.64279j] 10 [1.00000, 0.80902+0.58779j, 0.30902+0.95106j, -0.30902+0.95106j, -0.80902+0.58779j, -1.00000, -0.80902-0.58779j, -0.30902-0.95106j, 0.30902-0.95106j, 0.80902-0.58779j]
Ruby
Hopefully someone will fix the formatting <lang ruby>require 'complex'
for n in 2..10
printf "%2d ", n puts (0..n-1).map { |k| Complex.polar(1, 2 * Math::PI * k / n) }.join(" ")
end</lang> Output:
2 1.0+0.0i -1.0+1.22460635382238e-16i 3 1.0+0.0i -0.5+0.866025403784439i -0.5-0.866025403784438i 4 1.0+0.0i 6.12303176911189e-17+1.0i -1.0+1.22460635382238e-16i -1.83690953073357e-16-1.0i 5 1.0+0.0i 0.309016994374947+0.951056516295154i -0.809016994374947+0.587785252292473i -0.809016994374948-0.587785252292473i 0.309016994374947-0.951056516295154i 6 1.0+0.0i 0.5+0.866025403784439i -0.5+0.866025403784439i -1.0+1.22460635382238e-16i -0.5-0.866025403784438i 0.5-0.866025403784439i 7 1.0+0.0i 0.623489801858734+0.78183148246803i -0.222520933956314+0.974927912181824i -0.900968867902419+0.433883739117558i -0.900968867902419-0.433883739117558i -0.222520933956315-0.974927912181824i 0.623489801858733-0.78183148246803i 8 1.0+0.0i 0.707106781186548+0.707106781186547i 6.12303176911189e-17+1.0i -0.707106781186547+0.707106781186548i -1.0+1.22460635382238e-16i -0.707106781186548-0.707106781186547i -1.83690953073357e-16-1.0i 0.707106781186547-0.707106781186548i 9 1.0+0.0i 0.766044443118978+0.642787609686539i 0.17364817766693+0.984807753012208i -0.5+0.866025403784439i -0.939692620785908+0.342020143325669i -0.939692620785908-0.342020143325669i -0.5-0.866025403784438i 0.17364817766693-0.984807753012208i 0.766044443118978-0.64278760968654i 10 1.0+0.0i 0.809016994374947+0.587785252292473i 0.309016994374947+0.951056516295154i -0.309016994374947+0.951056516295154i -0.809016994374947+0.587785252292473i -1.0+1.22460635382238e-16i -0.809016994374948-0.587785252292473i -0.309016994374948-0.951056516295154i 0.309016994374947-0.951056516295154i 0.809016994374947-0.587785252292473i
Seed7
$ include "seed7_05.s7i"; include "float.s7i"; include "complex.s7i"; const proc: main is func local var integer: n is 0; var integer: k is 0; begin for n range 2 to 10 do write(n lpad 2 <& ": "); for k range 0 to pred(n) do write(polar(1.0, 2.0 * PI * flt(k) / flt(n)) digits 4 lpad 15 <& " "); end for; writeln; end for; end func;
Output:
2: 1.0000+0.0000i -1.0000+0.0000i 3: 1.0000+0.0000i -0.5000+0.8660i -0.5000-0.8660i 4: 1.0000+0.0000i 0.0000+1.0000i -1.0000+0.0000i 0.0000-1.0000i 5: 1.0000+0.0000i 0.3090+0.9511i -0.8090+0.5878i -0.8090-0.5878i 0.3090-0.9511i 6: 1.0000+0.0000i 0.5000+0.8660i -0.5000+0.8660i -1.0000+0.0000i -0.5000-0.8660i 0.5000-0.8660i 7: 1.0000+0.0000i 0.6235+0.7818i -0.2225+0.9749i -0.9010+0.4339i -0.9010-0.4339i -0.2225-0.9749i 0.6235-0.7818i 8: 1.0000+0.0000i 0.7071+0.7071i 0.0000+1.0000i -0.7071+0.7071i -1.0000+0.0000i -0.7071-0.7071i 0.0000-1.0000i 0.7071-0.7071i 9: 1.0000+0.0000i 0.7660+0.6428i 0.1736+0.9848i -0.5000+0.8660i -0.9397+0.3420i -0.9397-0.3420i -0.5000-0.8660i 0.1736-0.9848i 0.7660-0.6428i 10: 1.0000+0.0000i 0.8090+0.5878i 0.3090+0.9511i -0.3090+0.9511i -0.8090+0.5878i -1.0000+0.0000i -0.8090-0.5878i -0.3090-0.9511i 0.3090-0.9511i 0.8090-0.5878i
Scheme
<lang scheme>(define pi (* 4 (atan 1)))
(do ((n 2 (+ n 1)))
((> n 10)) (display n) (do ((k 0 (+ k 1))) ((>= k n)) (display " ") (display (make-polar 1 (* 2 pi (/ k n))))) (newline))</lang>
Tcl
<lang Tcl>package require Tcl 8.5 namespace import tcl::mathfunc::*
set pi 3.14159265 for {set n 2} {$n <= 10} {incr n} {
set angle 0.0 set row $n: for {set i 1} {$i <= $n} {incr i} { lappend row [format %5.4f%+5.4fi [cos $angle] [sin $angle]] set angle [expr {$angle + 2*$pi/$n}] } puts $row
}</lang>