Roots of unity

From Rosetta Code

Jump to: navigation, search

Programming Task
This is a programming task. It lays out a problem which Rosetta Code users are encouraged to solve, using languages they know.

Code examples should be formatted along the lines of one of the existing prototypes.
The purpose of this task is to explore working with complex numbers. Given n, find the n-th roots of unity.

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

This example is incorrect. It should print all n roots, either individually or by explicitly labeling conjugate pairs (x ± yi). Please fix the code and remove this message. 


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))
Retrieved from "http://rosettacode.org/wiki/Roots_of_unity"
Personal tools