✍ diale.org

Webpage of Tiago Charters de Azevedo

Início/Start Arquivo/Archive AdNauseum Notas/Notes RSS


Gráfico da função gamma

... para impressão 3d

2018/02/08-12:33:48

Tudo começou por um tuite sobre o gráfico de 1909 da extensão ao plano complexo da função gamma

Uns minutos depois fiz o upload para o Github do STL

e 10h depois alguém do outro do Atlântico tinha uma o modelo impresso.

Giro!

Etiquetas/Tags: 3d print, reprap, matemática

Hexy

a light-seeking robot – G. Draper (British)

2018/02/07-14:59:54

Note: great idea for an Arduino robot ;)

Etiquetas/Tags: hexy, robot, analog, 1965, diy, arduino

Some notes on a V-Bot

... math and physics

2018/01/11-13:49:57

A V-Bot is a simple drawing machine based on two-center bipolar coordinates. Here's some web references:

The layout of a V-Bot is quite simple, two motors at the upper vertices at a distance d, with two strings attached with lengths l1 and l2. The drawing pen rests at (x,y).

Problems to solve

The typical math problems associated with robots are two: direct and inverse kinematics. In this particular case the inverse kinematics is easy to obtain, given the pen position obtain the lengths of the strings:

latex2png equation

I'm using GNU/Octave, so here's some auxiliary functions:

function retval=l1c(x,y,d)
  retval= sqrt(y.^2+(x+d/2).^2);
end
function retval=l2c(x,y,d)
  retval= sqrt(y.^2+(x-d/2).^2);
end

These equations define a map from the cartesian coordinates (x,y) to (l1,l2).

For the direct kinematics one has, with some simple algebra, latex2png equation

Drawing resolution

One can tackle the problem of finding the resolution in two ways. Using the Jacobian and the arc length of a curve.

The Jacobian

The Jacobian latex2png equation defines the ratio of a unit area in both coordinates. That is, it controls what happens (direct kinematics) when a unit area in (l1,l2) coordinates is transformed in to (x,y) coordinates, it expands, contracts or stays the same. Values for J greater than one gives an area expansion, lower than one the area is reduced. The Jacobian controls the resolution of the drawing region.

What happens when the Jacobian is zero? In simple terms it simply means that the transformation breaks, i.e. one can not use the transformation between (x,y) and (l1,l2) to control the robot, and so at the points or at the lines where Jacobian is zero the robot can not be controlled (you can guess by simple inspection, with no math, what these lines/points are ;) ).

The Jacobian latex2png equation is given by latex2png equation

The next image shows the value of this Jacobian as a function of (x,y)

The code for the Jacobian J(x,y,l1,l2) is given by

function retval=jac(x,y,d)
  retval=l1c(x,y,d).*l2c(x,y,d)./(d*y);
end

For the inverse kinematics the Jacobian is the reciprocal of J(x,y,l1,l2), that is 1/J(x,y,l1,l2). In this case one gets

As expected (see pics) the problematic lines are l1+l2=d or simply the line y=0 for any x.

The previous plots shows the "good" areas for drawing, somewhere below the y=0 line (10cm?). Will see the same result when considering the tension on the strings (see below).

The arc length of a curve

The the arc length of a curve can be used to determine the "good" plotting region for the robot. This is done by taking the first order Taylor expansion of the direct and inverse kinematics equations. This however requires some prior considerations.

Because the relation between (x,y) and (l1,l2) is non-linear any variation (x,y) will produce a non-linear variation in (l1,l2) this variation is path dependent, the propagation of the variations depends if, for example, the pen goes right and then up or by the hypotenuse to the same point on the drawing board.

latex2png equation and latex2png equation

Here goes:

function retval=lenghtl(x,y,dx,dy,d)
  dl1=(y.*dy+(x+d/2).*dx)./l1c(x,y,d);
  dl2=(y.*dy+(x-d/2).*dx)./l2c(x,y,d);
  retval=sqrt((dl1.^2+dl2.^2)/(dx.^2+dy.^2));
end
function retval=lenghtxy(x,y,dl1,dl2,d)
  dx=l1c(x,y,d)./d.*dl1 - l2c(x,y,d)./d.*dl2;
  dy=l1c(x,y,d).*(1-2*x/d)./(2*y).*dl1 + l2c(x,y,d).*(1+2*x/d)./(2*y).*dl2;
  retval=sqrt((dx.^2+dy.^2)/(dl1.^2+dl2.^2));
end

The next plot shows the ratio of unit lengths in the direct kinematics plotted in the (x,y) plane by taking the average on the 4 movements of the pen (right,0), (left,0), (0,up) and (0,down)1.

The next plot takes the average length on 4 movements of the pen (left,up), (right,up), (left,down) and (right,down)2.

This last plot is similar to the final plot obtained by Bill Rasmussen. My point is that in this way we only account for the variation of the length of a straight segment in this coordinate transformation. Obviously the length changes after the transformation due to the non-linearity of the transformation and reflects the loss of resolution but the right way to determine the resolution is using the Jacobian.

The ends of the control lines in the above picture seem to be further away from the plotting surface than V plotters commonly seen on the internet.

That's just because the important tool to measure resolution is the Jacobian, even thought many of the V-Bot builders do not know about it ;)

What about the tension on the strings?

One should also take into account the tension on each string. With some trigonometry one gets

latex2png equation

There's only one singular configuration that one should take note, the case latex2png equation which yields latex2png equation

function retval=tension(x,y,d)
  m=1;
  g=1;
  l1=sqrt(y.^2+(d/2-x).^2);
  l2=sqrt(y.^2+(d/2+x).^2);       
  cosa1=(d/2+x)./l1;
  cosa2=(d/2-x)./l2;
  sina1=y./l1;
  sina2=y./l2;
  
  T1=m*g*cosa2./(cosa1.*sina2+cosa2.*sina1);
  T2=m*g*cosa1./(cosa1.*sina2+cosa2.*sina1);
  retval=sqrt(T1.^2+T2.^2);
end

This is shown in the following plot:

Also the cases of null tension on one of the strings latex2png equation, latex2png equation or latex2png equation, latex2png equation yields, respectivly, latex2png equation and latex2png equation or latex2png equation and latex2png equation.

1. Vertical displacements: (right,0), (left,0), (0,up) and (0,down)

2. Oblique displacements: (left,up), (right,up), (left,down) and (right,down)

Etiquetas/Tags: v-bot, dc motor, math, physics, polargraph

Erros científicos nos meios de counicação social em Portugal

... uma lista não exaustiva

2017/10/15-12:33:45

Hortas inteligentes e aquecedores à luz das velas

Detetar ondas gravitacionais a partir da Terra ou do espaço

A construção de uma plotter reciclando os motores passa-a-passo de umas drives de CDROM não é um projecto original mas coloca alguns desafios de montagem, calibração e programação.

A minha abordagem não difere de muitas outras que se podem encontrar pela net fora, apenas difere no objectivo. Desenhar numa folha de papel, nem que seja de dimensões reduzidas (4cm por 4cm), soluções de equações diferenciais ordinárias.

Assim para além de ter de controlar os motores passo-a-passo é necessário construir a rotina de integração numérica. Implementei um Runge-Kutta de ordem 4 (há quantos exemplos na net) apenas tive de o ajustar às minhas necessidades.

A construção ficou assim:

Está disponível um vídeo.

As funções que implementam a resolução numérica são duas.

// Runge-Kutta

float rk(float t, float dt, float x[]){
  float xaux[]={0,0,0};
  int i=0,j=0;

    f(t,x,k1);
    for(i=0;i<m;i++){
        xaux[i]=x[i]+dt/2*k1[i];}

    f(t+dt/2,xaux,k2);
    for(i=0;i<m;i++){
      xaux[i]=x[i]+dt/2*k2[i];}

    f(t+dt/2,xaux,k3);
    for(i=0;i<m;i++){
      xaux[i]=x[i]+dt/2*k3[i];}
    
    f(t+dt,xaux,k4);

    for(i=0;i<m;i++){
      x[i]=x[i]+dt/6.0*(k1[i]+2.0*k2[i]+2*k3[i]+k4[i]);}

    t=t+dt;}

void f(float t, float x[],float aux[]){
  aux[0]=c1*(x[1]-x[0]);
  aux[1]=x[0]*(c2-x[2]);
  aux[2]=x[0]*x[1]-c3*x[2];}

E o código completo é este

#include <Servo.h>
#include "AFMotor.h"
#define LINE_BUFFER_LENGTH 512

/* -------------------------------------------------- */
// Plotter parameters
// SINGLE, DOUBLE. INTERLEAVE or MICROSTEP.
char STEP=MICROSTEP;
const int stepsPerRevolution=20; 
AF_Stepper PStepperY(stepsPerRevolution,1);            
AF_Stepper PStepperX(stepsPerRevolution,2);  

float StepsPerMillimeterX=100.0;
float StepsPerMillimeterY=100.0;
const float pi=3.14159265358979;

const int penZUp=80;
const int penZDown=66;
const int penServoPin =10; 
Servo penServo;  
/* -------------------------------------------------- */

/* -------------------------------------------------- */
// Drawing robot limits, in mm
float Xmin=0;
float Xmax=40;
float Ymin=0;
float Ymax=40;
float Zmin=0;
float Zmax=1;

float Xpos=Xmin;
float Ypos=Ymin;
float Zpos=Zmax; 

float StepInc=1;
int StepDelay=0;
int LineDelay=0;
int penDelay=5;
/* -------------------------------------------------- */

/* -------------------------------------------------- */
// begin ode parameters
const float c1=10.0; 
const float c2=28.0;
const float c3=8.0/3.0;
const int m=3; // ode order
const float dt=0.01; //integration step
  
float x[m]={1, 1, 1};    // init conditions
float xaux[m]={0, 0, 0};    // aux
float k1[m], k2[m], k3[m], k4[m]; // arrays for the Runge-Kutta intermediate points
/* -------------------------------------------------- */

void setup() {
  //  Serial.begin(9600);

  PStepperX.setSpeed(5);
  PStepperY.setSpeed(5);
  
  penServo.attach(penServoPin);
  penServo.write(penZUp);
  delay(100);}

void loop(){
  int i=0;
  int n=10000;
  float t=0;

  penUp();
  delay(1000);

  /* -------------------------------------------------- */
  // Do some jogging...
  /* -------------------------------------------------- */
  move(Xmax,0);
  move(Xmax,Ymax);
  move(0,Ymax);
  move(0,0);
  /* -------------------------------------------------- */

  
  /* -------------------------------------------------- */
  // Draw the square limits
  /* -------------------------------------------------- */
  drawLine(Xmax,0);
  drawLine(Xmax,Ymax);
  drawLine(0,Ymax);
  drawLine(0,0);
  /* -------------------------------------------------- */
  
  penUp();
  delay(2000);

  /* -------------------------------------------------- */
  // Move to initial position
  /* -------------------------------------------------- */
  move(x[1],x[2]);

  
  // use random initials conditions
  // or close to the one given {1,1,1}
  for(i=0;i<=n;i++){
    rk(t,dt, x);
    drawLine(Xmax/2+x[1]*2.0/3.0,x[2]*4.0/5.0);
  }

  
    penUp();
    move(Xmax,Ymax);
    delay(100000);}

void penUp(){ 
   penServo.write(penZUp); 
   delay(penDelay);}
 
 void penDown(){ 
   penServo.write(penZDown); 
   delay(penDelay);}
 
void drawLine(float x1, float y1){
  if (x1 >= Xmax){x1=Xmax;}
  if (x1 <= Xmin){x1=Xmin;}
  if (y1 >= Ymax){y1=Ymax;}
  if (y1 <= Ymin){y1=Ymin;}
  
  x1=(int)(x1*StepsPerMillimeterX);     //  Convert coordinates to steps
  y1=(int)(y1*StepsPerMillimeterY);
  float x0=Xpos;
  float y0=Ypos;
  
  long dx=abs(x1-x0);     //  Let's find out the change for the coordinates
  long dy=abs(y1-y0);
  int sx=x0<x1 ? StepInc : -StepInc;
  int sy=y0<y1 ? StepInc : -StepInc;
  
  long i;
  long over=0;
  penDown();
  
  if (dx>dy){
    for (i=0; i<dx; ++i) {
      PStepperX.onestep(sx,STEP);
      over+=dy;
      if (over>=dx) {
        over-=dx;
        PStepperY.onestep(sy,STEP);}
      delay(StepDelay);}}
  else{
    for (i=0; i<dy; ++i){
      PStepperY.onestep(sy,STEP);
      over+=dx;
      if (over>=dy){
        over-=dy;
        PStepperX.onestep(sx,STEP);}
      delay(StepDelay);}}
  
  Xpos=x1;
  Ypos=y1;}


void move(float x1, float y1){
  penUp();
  
  if (x1 >= Xmax){
    x1=Xmax;}
  if (x1 <= Xmin){
    x1=Xmin;}
  if (y1 >= Ymax){
    y1=Ymax;}
  if (y1 <= Ymin){
    y1=Ymin;}
  
  x1=(int)(x1*StepsPerMillimeterX);     //  Convert coordinates to steps
  y1=(int)(y1*StepsPerMillimeterY);
  float x0=Xpos;
  float y0=Ypos;
  
  long dx=abs(x1-x0);     //  Let's find out the change for the coordinates
  long dy=abs(y1-y0);
  int sx=x0<x1 ? StepInc : -StepInc;
  int sy=y0<y1 ? StepInc : -StepInc;
  
  long i;
  long over=0;
  
  if (dx>dy){
    for (i=0; i<dx; ++i) {
      PStepperX.onestep(sx,STEP);
      over+=dy;
      if (over>=dx) {
        over-=dx;
        PStepperY.onestep(sy,STEP);}
      delay(StepDelay);}}
  else{
    for (i=0; i<dy; ++i){
      PStepperY.onestep(sy,STEP);
      over+=dx;
      if (over>=dy){
        over-=dy;
        PStepperX.onestep(sx,STEP);}
      delay(StepDelay);}}
  
  Xpos=x1;
  Ypos=y1;}


// Runge-Kutta

float rk(float t, float dt, float x[]){
  float xaux[]={0,0,0};
  
  int i=0,j=0;
  
    
    f(t,x,k1);
    for(i=0;i<m;i++){
        xaux[i]=x[i]+dt/2*k1[i];}

    f(t+dt/2,xaux,k2);
    for(i=0;i<m;i++){
      xaux[i]=x[i]+dt/2*k2[i];}

    f(t+dt/2,xaux,k3);
    for(i=0;i<m;i++){
      xaux[i]=x[i]+dt/2*k3[i];}
    
    f(t+dt,xaux,k4);

    for(i=0;i<m;i++){
      x[i]=x[i]+dt/6.0*(k1[i]+2.0*k2[i]+2*k3[i]+k4[i]);}

    t=t+dt;

  }

void f(float t, float x[],float aux[]){
  aux[0]=c1*(x[1]-x[0]);
  aux[1]=x[0]*(c2-x[2]);
  aux[2]=x[0]*x[1]-c3*x[2];}

Parte das peças usadas foram impressas e modeladas usando o OpenScad:

include <../utils/polyholes.scad>
$fn=64;
h=15;


module pin(r=2.8){
    difference(){
        union(){
            cylinder(h=h,r=5,center=true);
            translate([0,0,1]){
                cylinder(h,r=r,center=true);}}
        poly_cylinder(h=2*h,r=2,center=true);}}

module print(){
    n=3;
    translate([0,0,h/2]){
        for(i=[0:n]){
            translate([0,12*i,0]){
                pin(2.8);}
            translate([12,12*i,0]){
                pin(7.1/2);}}}

    translate([20,0,0]){
        rotate([0,0,90])
        fixing();}

        translate([20,11*2.5,0]){
        rotate([0,0,90])
        fixing();}}



module fixing(l=25,rc=.1*22){
    difference(){
        minkowski(){
            cube([l-rc,l-rc,l-rc],center=true);
            sphere(rc);}
        
        translate([0,0,-l/2+.01]){
            cube([2*l,2*l,l],center=true);}
        
        translate([0,l/2,l]){
            cube([2*l,l,2*l],center=true);}
        
        translate([0,-l/1.4,l/1.4]){
            rotate([45,0,0])
            cube([1.5*l,l,l],center=true);}
        
        translate([0,-l/4,l-4.5]){
            cylinder(l,r=4,center=true,$fn=6);
            poly_cylinder(2*l,r=2,center=true);}
        
        translate([l/3,0,l/4]){
            rotate([90,0,0]){
                translate([0,0,l-4.5]){
                    cylinder(l,r=4,center=true,$fn=6);
                    poly_cylinder(2*l,r=2,center=true);}}}

        translate([-l/3,0,l/4]){
            rotate([90,0,0]){
                translate([0,0,l-4.5]){
                    cylinder(l,r=4,center=true,$fn=6);
                    poly_cylinder(2*l,r=2,center=true);}}}}}    


module otherfixing(l=25,rc=.1*22){
    difference(){
        minkowski(){
            cube([l-rc,l-rc,l-rc],center=true);
            sphere(rc);}
        
        translate([0,0,-l/2+.01]){
            cube([2*l,2*l,l],center=true);}
        
        translate([0,l/2,l]){
            cube([2*l,l,2*l],center=true);}
        
        translate([0,-l/1.4,l/1.4]){
            rotate([45,0,0])
            cube([1.5*l,l,l],center=true);}
        
        translate([0,-l/4,l-4.5]){
            cylinder(l,r=4.5,center=true,$fn=64);
            poly_cylinder(2*l,r=2,center=true);}
        
        translate([l/3,0,l/4]){
            rotate([90,0,0]){
                translate([0,0,l-4.5]){
                    cylinder(l,r=4.5,center=true,$fn=66);
                    poly_cylinder(2*l,r=2,center=true);}}}

        translate([-l/3,0,l/4]){
            rotate([90,0,0]){
                translate([0,0,l-4.5]){
                    cylinder(l,r=4.5,center=true,$fn=64);
                    poly_cylinder(2*l,r=2,center=true);}}}}}    



a=25.4*4/3;
b=25.4*2;

module penbase(){
    
    difference(){
        cube([a,b,2.5],center=true);
        translate([0,12,10/2]){
            rotate([0,0,90]){
                poly_cylinder(h=b,r=2.1,center=true);}}}


    //Servo holder
    difference(){
        translate([a/2-5/2,-b/2+2.5,15/2-2.5/2]){
            cube([5,5,15],center=true);}
        
        translate([a/2-1.25/2,-b/2+5-2,2+5]){
            rotate([0,90,0]){
                cylinder(10,r=.8,center=true);}}}

    //Servo holder
    difference(){
        translate([a/2-5/2,-b/2+5+2.5+23.5,15/2-2.5/2]){
            cube([5,5,15],center=true);}
        
        translate([a/2-1.25/2,-b/2+21+7+2.5,2+5]){
            rotate([0,90,0]){
                cylinder(h=10,r=.8,center=true);}}}
    
    difference(){
        union(){
            translate([-a/2+6,b/2-2.5,10/2-2.5/2]){
                cube([12,5,10],center=true);}
            
            translate([-a/2+6,-b/2+2.5,10/2-2.5/2]){
                cube([12,5,10],center=true);}}
        
        
        translate([-a/2+6-3,2.5,10/2]){
            rotate([0,90,90]){
                poly_cylinder(h=b,r=3.1/2,center=true);}}
        
        translate([-a/2+6+3,2.5,10/2]){
            rotate([0,90,90]){
                poly_cylinder(h=b,r=3.1/2,center=true);}}


    }}

module penholder(){
    difference(){
        translate([-a/2+6,0,22/2+2.5/2+1]){
            cube([13,12,22],center=true);}

        translate([-a/2+6-3,2.5,10/2]){
            rotate([0,90,90]){
                poly_cylinder(h=b,r=3.3/2,center=true);}}

        translate([-a/2+6+3,2.5,10/2]){
            rotate([0,90,90]){
                poly_cylinder(h=b,r=3.3/2,center=true);}}
        
        translate([-a/2+6,2.5,13]){
            rotate([0,90,90]){
                poly_cylinder(h=b,r=4.25,center=true);}}
        
        translate([-a/2+6,2.5,25]){
            rotate([0,0,90]){
                cube([20,2,20],center=true);}}
        
        translate([-a/2+6,0,21]){
            rotate([0,90,0]){
                poly_cylinder(h=b,r=2,center=true);}}
        
        
    }}

module otherpin(h=7,r=2){
    difference(){

            cylinder(h=h,r=5,center=true);      
        poly_cylinder(h=2*h,r=2,center=true);}}


module foot(h=12,r=2){
    difference(){
        cylinder(h=h,r1=5,r2=5,center=true);
        translate([0,0,-h+5]){
            poly_cylinder(h,r=6.25/2,center=true,$fn=6);}
        poly_cylinder(h=1.1*h,r=2.1,center=true);}}

module penpin(h=7,r=2){
    difference(){
        hull(){
            cylinder(h=h,r=5,center=true);
            translate([0,5,0])
            cylinder(h=h,r=5,center=true);}

        
        poly_cylinder(h=2*h,r=r,center=true);}}

 penbase();

Etiquetas/Tags: dvd, plotter, maker, diy, arduino

Isto é parte da dificuldade, apenas o início.

Senso comum e a sua sua falta dele

2017/09/20-12:06:51

Etiquetas/Tags: filosofia, física, senso comum

Yet another Platonic solid set

... circumscribed in a 25mm radius sphere

2017/05/25-12:35:21

Get it here: https://www.youmagine.com/designs/yet-another-platonic-solid-set

Print it, share and learn!

Etiquetas/Tags: 3dprint, 3dprinting, Platonic, regular, solid, solids, polyhedron

So you want to print a Truncated Icosahedron ?

Buckyball

2017/05/25-12:29:36

Get it here: https://www.youmagine.com/designs/truncated-icosahedron-buckyball

Print it, share and learn!

Etiquetas/Tags: buckyball, Truncated Icosahedron, Truncated, Icosahedron, 3dprint, 3dprinting

Paradroid: how I did it

Steps to build a 2.5D Paradroid for 3D printing

2017/05/03-14:11:33

It all started when I found some drawings of the character of the infamous computer game from the 80's.

The next step was obvious to me, just build a raster graphics of the beast using GNU/Octave. As is the original image needed some modification on order to be 3d-printed.

A=[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1;
1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1;    
1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1;
1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1;
1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;    
1 1 0 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1;
1 0 0 0 0 0 1 1 1 0 0 0 0 0 1 1 1 1 0 0 0;
0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0;
0 0 1 1 1 0 0 1 0 0 1 1 1 0 0 1 1 1 0 0 0;
0 0 1 1 1 0 0 1 0 0 1 1 1 0 0 1 1 1 0 0 0;
0 0 1 1 1 0 0 1 0 0 1 1 1 0 0 1 1 1 0 0 0;
0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 1 1 0 0 0;
1 0 0 0 0 0 1 1 1 0 0 0 0 0 1 1 1 1 0 0 0;    
1 1 0 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1;
1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;
1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1;    
1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1];
imshow(A)

print -deps -color paradroid.eps
print -dpng -color paradroid.png

So added some extra squares to get the job done.

Then used the Inkscape-to-OpenSCAD plugin and got the final 3d model, after some editing on the output code.

And it prints nicely with a resolution of 0.2mm and a 20% infill.

You can get it here: https://www.youmagine.com/designs/paradroid-ornament-key-chain

Share and enjoy!

Etiquetas/Tags: paradroid, 3d, 3dprinting, 3dprint, game, math

Números pseudo-aleatórios em LISP

Implementação em LISP de um gerador de números pseudo-aleatórios

2017/04/19-17:34:32

Uma das formas de gerar números com uma distribuição uniforme num intervalo é usando o método de congruências lineares1. Vejamos então como implementar um gerador de números pseudo-aleatórios que gera números inteiros entre 0 e m com uma distribuição uniforme latex2png equation.

A ideia é gerar uma sucessão de números com a forma latex2png equation onde2

latex2png equation

Ao valor de latex2png equation dá-se o nome de semente porque é o primeiro termo da sucessão e é usado como input para o termo seguinte. Note-se que para um mesmo latex2png equation obtém-se sempre o mesmo latex2png equation.

A sucessão de números assim obtida tem período latex2png equation e o intervalo de variação dos valores da sucessão pode ser ajustado ao intervalo latex2png equation através de

latex2png equation

A maneira mais directa é escrever

(defun rand (a c m x)
  (mod (+ (* a x) c) m))
e obtém-se
> (rand 24298 99991 199017 0)

99991

Esta não é a maneira preferível de obter um número aleatório. Queremos obter um número simplesmente fazendo (rand), mas para isso, temos de ter uma maneira de guardar os sucessivos termos da sucessão de modo a que sirvam de sementes para os termos seguintes. Assim vamos definir a variável

(defvar *r-seed* 0)
que inicializa a sucessão de sementes (ou dos próprios números pseudo-aleatórios). Esta é uma variável global e por isso está enquadrada por *. Definimos um substituto para a função anterior como
(defun rand1 ()
  (let ((a 24298)
        (c 99991)
        (m 199017))
    (cond ((<= *r-seed*  199017)
           (setf *r-seed* (mod (+ (* a *r-seed*) c) m)))
          (t
           (setf *r-seed* 0)))))
onde agora os parâmetros a, c e m são locais e estão definidos dentro da função. Assim
> (rand1)

99991

A função rand1 é simples de perceber. Depois de inicializar as quantidades a, c e m verifica se *r-seed* é menor ou igual a 199017, se sim altera o valor de *r-seed* para o novo termo da sucessão através de

(setf *r-seed* (mod (+ (* a *r-seed*) c) m))
se não
(setf *r-seed* 0)

De modo a gerar números aleatórios entre latex2png equation define-se a função

(defun rand (&optional alpha beta)
  (let ((m 199017))
    (cond ((not (and alpha beta))
           (float (/ (rand1) m)))
          (t
           (+ (* (float (/ (rand1) m)) (- beta alpha)) alpha)))))
com dois argumentos opcionais alpha e beta; se não forem dados (rand) gera números com uma distribuição uniforme no intervalo latex2png equation:
> (rand)

0.55237997

> (rand 0 2)

1.1378727

Ref:

1. D. Knuth, TAOCP

2. Master Library da TI Programable 58/59

Etiquetas/Tags: LISP, números pseudo-aleatórios

Palavras chave/keywords: página pessoal, blog

Criado/Created: NaN

Última actualização/Last updated: 08-02-2018 [12:34]


GNU/Emacs

1999-2017 (ç) Tiago Charters de Azevedo

São permitidas cópias textuais parciais/integrais em qualquer meio com/sem alterações desde que se mantenha este aviso.

Verbatim copying and redistribution of this entire page are permitted provided this notice is preserved.