Roots of unity
From Rosetta Code
Programming Task
This is a programming task. It lays out a problem which Rosetta Code users are encouraged to solve, using languages they know.
Contents |
[edit] 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;
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
[edit] 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
[edit] BASIC
Works with: QuickBasic version 4.5
Translation of: Java
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
[edit] 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"); } }
[edit] C++
#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; } }
[edit] Common Lisp
(defun roots-of-unity (n)
(loop for i below n
collect (cis (* pi (/ (* 2 i) n)))))
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)).
[edit] D
Works with: D version 2.012
Works with: D version 1.028
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)) ; }
[edit] 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
[edit] Fortran
Works with: Fortran version 90 and later
PROGRAM Roots
COMPLEX :: root
INTEGER :: n
REAL :: angle, pi
pi = 4.0 * ATAN(1.0)
DO n = 2, 7
angle = 0.0
WRITE(*,"(I1,A)", ADVANCE="NO") n,": "
DO WHILE (angle < 6.283)
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
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
[edit] 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)]
[edit] 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)
[edit] 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 )
[edit] 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.
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 } }
[edit] 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
[edit] Perl
Works with: Perl version 5.8.8
Library: Math::Complex
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"; }
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
[edit] Python
Works with: Python version 2.5.1
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)) for nr in range(2,11): print nr, list(croots(nr))
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]
[edit] 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
[edit] 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))
Categories: Programming Tasks | Arithmetic operations | Ada | ALGOL 68 | ALGOL 68 examples needing attention | Examples needing attention | BASIC | C | C++ | Common Lisp | D | Forth | Fortran | Haskell | IDL | J | Java | OCaml | Perl | Math::Complex | Python | Seed7 | Scheme

