Re : Function vs list of instructions executed separately

["Ewan Delanoy" <ewan.delanoy@gmx.fr>]

>1) You almost certainly want to use my() and not local() 2) If this does not fix the problem, can you >send the whole script script to me, I'd like to be able to run all commands :-)>

 I changed all the local’s to my’s in my code, and the problem stays the same.
  So I append below my ugly, personal, unreadable code (it compiles fine though).
(in that new version, items number %26 and %27 should be equal and they are not).

Ewan

/************************************************/

eb_variables_array=['eb1,'eb2,'eb3,'eb4,'eb5,'eb6,'eb7,'eb8,'eb9,'eb10,'eb11,'eb12,'eb13,'eb14,'eb15,'eb16,'eb17,'eb18,'eb19,'eb20,'eb21,'eb22,'eb23,'eb24,'eb25,'eb26,'eb27,'eb28,'eb29,'eb30,'eb31,'eb32,'eb33,'eb34,'eb35,'eb36,'eb37,'eb38,'eb39,'eb40,'eb41,'eb42,'eb43,'eb44,'eb45,'eb46,'eb47,'eb48,'eb49,'eb50,'eb51,'eb52,'eb53,'eb54,'eb55,'eb56,'eb57,'eb58,'eb59,'eb60,'eb61,'eb62,'eb63,'eb64,'eb65,'eb66,'eb67,'eb68,'eb69,'eb70,'eb71,'eb72,'eb73,'eb74,'eb75,'eb76,'eb77,'eb78,'eb79,'eb80,'eb81,'eb82,'eb83,'eb84,'eb85,'eb86,'eb87,'eb88,'eb89,'eb90,'eb91,'eb92,'eb93,'eb94,'eb95,'eb96,'eb97,'eb98,'eb99,'eb100]

ryan_valuation(expr,l_var)=
{
   my(n,k,answer);
   n=length(l_var);
   answer=0;
   k=1;
   while((answer<1)*(k<=n),if(polcoeff(expr,1,l_var[k])!=0,answer=k;);k++;);
   return(answer);
}

ryan_degree(expr,l_var)=
{
   my(n,k,answer);
   n=length(l_var);
   answer=0;
   k=n;
   while((answer<1)*(k>0),if(polcoeff(expr,1,l_var[k])!=0,answer=k;);k--;);
   return(answer);

}


first_nonzero_ryan_degree(arr,l_var)=
{
  my(n,k,answer);
   n=length(arr);
   answer=0;
   k=1;
   while((answer<1)*(k<=n),if(ryan_degree(arr[k],l_var)>0,answer=k;);k++;);
   return(answer);
}


compute_kernel(l_expr,l_var)=
{
   my(variable_l,watcher,answer,part1,part2,expr,var1,form1,\
   new_l_var,unchanged_indices,unchanged_variables,n,k0,k1);
   part1=[];
   part2=[];
   variable_l=l_expr;
   watcher=1;
   n=length(l_expr);
   while(watcher,
      k0=first_nonzero_ryan_degree(variable_l,l_var);
      if(k0==0,watcher=0,
         expr=variable_l[k0];
         k1=ryan_degree(expr,l_var);
         var1=l_var[k1];
         form1=make_zero(expr,var1);
         part1=concat(part1,[var1]);
         part2=concat(part2,[form1]);
         variable_l=subst(variable_l,var1,form1);
      );
   );
   new_l_var=big_subst(l_var,part1,part2);
   unchanged_indices=List([]);
   unchanged_variables=List([]);
   annulators=List([]);
   for(k=1,length(l_var),\
   if(l_var[k]==new_l_var[k],\
   listput(unchanged_indices,k);\
   listput(unchanged_variables,l_var[k]);,\
   listput(annulators,l_var[k]-new_l_var[k]);\
   ););
   dual_version=vector(length(unchanged_variables),k,\
   sum(j=1,length(l_var),polcoeff(new_l_var[j],1,unchanged_variables[k])*l_var[j]));
   return([new_l_var,Vec(annulators),dual_version,Vec(unchanged_indices),Vec(unchanged_variables)]);

}



extract_basis(l_expr,l_var)=
{
    my(rewritings,counter,current_var,current_renamed_var,\
    current_expr,writings,n,k,d,cj,bk,ck,rk);
    n=length(l_expr);
    rewritings=vector(n,k,[0,l_expr[k]]);
    counter=0;
    writings=List([]);
    for(k=1,n,\
    current_expr=rewritings[k][2];\
    d=ryan_degree(current_expr,l_var);\
    if(d>0,\
    current_var=l_var[d];\
    counter++;\
    current_renamed_var=eb_variables_array[counter];\
    listput(writings,counter);
    bk=rewritings[k][1];\
    ck=polcoeff(current_expr,1,current_var);\
    rk=current_expr-ck*current_var;
    for(j=k,n,\
    cj=polcoeff(rewritings[j][2],1,current_var);\
    rewritings[j][1]-=(cj/ck)*bk;\
    rewritings[j][1]+=(cj/ck)*current_renamed_var;\
    rewritings[j][2]-=cj*current_var;\
    rewritings[j][2]-=(cj/ck)*rk;\
    );););
    return([Vec(writings),rewritings]);
}

compute_rank(l_expr,l_var)=length(extract_basis(l_expr,l_var)[1])



big_concat(l)={\
my(answer,k);
answer=l[1];
for(k=2,length(l),answer=concat(answer,l[k]));
return(answer);
}


next_uple(bounds,uple)={
 my(answer,k,n,watcher,largest_index_for_subcritical_number);
 n=length(uple);
 largest_index_for_subcritical_number=0;
 watcher=1;
 k=n;
 while(watcher,\
 if(uple[k]<bounds[k],\
 watcher=0;largest_index_for_subcritical_number=k;,\
 k--;watcher=(k>0););\
 );
 if(largest_index_for_subcritical_number==0,return(vector(n,j,1)););
 answer=uple;
 answer[k]+=1;
 for(j=k+1,n,answer[j]=1;);
 return(answer);
}
 


cartesian_big_product(list_of_lists)={
   my(d,k,sizes,main_size,indices,uple,answer);
   d=length(list_of_lists);
   sizes=vector(d,k,length(list_of_lists[k]));
   main_size=prod(k=1,d,sizes[k]);
   answer=List([]);
   indices=vector(d,k,1);
   for(j=1,main_size,\
   uple=vector(d,k,list_of_lists[k][indices[k]]);
   listput(answer,uple);
   indices=next_uple(sizes,indices);
   );
   return(answer);
}

cartesian_product(l1,l2)=cartesian_big_product([l1,l2])


scalar_product_according_to_list(vec1,vec2,list_of_vars)=\
sum(k=1,length(list_of_vars),\
polcoeff(vec1,1,list_of_vars[k])*polcoeff(vec2,1,list_of_vars[k]))


test_individual_for_extremality(old_extr_constraints,new_vector,list_of_vars)=
{
    my(k,r,accu);
    accu=List([]);
    for(k=1,length(old_extr_constraints),\
    if(scalar_product_according_to_list(old_extr_constraints[k],\
    new_vector,list_of_vars)==0,\
    listput(accu,old_extr_constraints[k]);)\
    );
    r=compute_rank(accu,list_of_vars);
    return(r==length(list_of_vars)-1);
}

test_individuals_for_extremality(old_extr_constraints,new_vectors,list_of_vars)=
{
    my(k,r,accu);
    accu=List([]);
    for(k=1,length(new_vectors),\
    if(test_individual_for_extremality(old_extr_constraints,new_vectors[k],list_of_vars),\
    listput(accu,new_vectors[k]);)\
    );
    return(Vec(accu));
}

untested_new_vectors(extr_constraints,extr_vectors,new_cstr,list_of_vars)=
{
   my(k,scal_prods,accu_neg_indices,accu_zer_indices,accu_pos_indices,\
   temp1,part1,part2,part3,whole);
   scal_prods=vector(length(extr_vectors),k,\
   simplify(scalar_product_according_to_list(new_cstr,extr_vectors[k],list_of_vars))\
   );
   accu_neg_indices=List([]);
   accu_zer_indices=List([]);
   accu_pos_indices=List([]);
   for(k=1,length(scal_prods),\
   if(scal_prods[k]<0,listput(accu_neg_indices,k););\
   if(scal_prods[k]==0,listput(accu_zer_indices,k););\
   if(scal_prods[k]>0,listput(accu_pos_indices,k););\
   );
   part1=vector(length(accu_zer_indices),k,extr_vectors[accu_zer_indices[k]]);
   part2=vector(length(accu_pos_indices),k,extr_vectors[accu_pos_indices[k]]);
   temp1=cartesian_product(accu_pos_indices,accu_neg_indices);
   part3=vector(length(temp1),k,\
   (-1)*scal_prods[temp1[k][2]]*(extr_vectors[temp1[k][1]])+\
   scal_prods[temp1[k][1]]*(extr_vectors[temp1[k][2]])\
   );
   whole=big_concat([part1,part2,part3]);
   return(whole);
}



insert_constraint0(extr_constraints,extr_vectors,new_cstr,list_of_vars)=
{
   my(temp_vecs,temp_cstr,new_vecs,new_cstr);
   temp_vecs=untested_new_vectors(extr_constraints,extr_vectors,new_cstr,list_of_vars);
   temp_cstr=concat(extr_constraints,new_cstr);
   new_vecs=test_individuals_for_extremality(temp_cstr,temp_vecs,list_of_vars);
   new_cstr=test_individuals_for_extremality(new_vecs,temp_cstr,list_of_vars);
   return([new_cstr,new_vecs]);
}


extr_constraints=[i, j, k, u]
extr_constraints=[i, j, k, u]
extr_vectors=[i, j, k, u]
new_cstr=4*u-i
list_of_vars=[i,j,k,u]

temp_vecs=untested_new_vectors(extr_constraints,extr_vectors,new_cstr,list_of_vars)
temp_cstr=concat(extr_constraints,new_cstr)
new_vecs=test_individuals_for_extremality(temp_cstr,temp_vecs,list_of_vars)
new_cstr=test_individuals_for_extremality(new_vecs,temp_cstr,list_of_vars)
[new_cstr,new_vecs]

insert_constraint0(extr_constraints,extr_vectors,new_cstr,list_of_vars)