product except zero constraint

[Nartoo Meon]

Greetings!
How to implemet a constraint like V #= prod(DVars) but with only non-zero values from  DVars? A "prod_except_zero()" constraint.

I want to use it at a runtime, so it seems difficult to me. My previous tries always fail. Sorry I'm not so "professional".
Perhaps, it possible to use member()-like predicate and somehow save non-zero values while it update the target value.

Or use a cl_facts() but I don't know how.

[Hakan Kjellerstrand]

Hi Nartoo.

One idea is to use cond/3 to replace each 0 with 1 in the list and then using the standard prod/1 constraint.

For a list of only 0s this will yield 1, which I assume is correct.

"""

import cp.

main ?=>
  X = new_list(3),
  X :: 0..4,

  % At least one 0
  count(0,X,#>,0),

  Prod #= prod([cond(X[I]#=0,1,X[I]) : I in 1..X.len]), 

  % product_except_0(X,Prod),

  solve(X ++ [Prod]),
  println(X=Prod),
  fail,
  nl.
main => true.

% As a predicate
product_except_0(X,Prod) =>
  Prod #= prod([cond( X[I]#=0,1,X[I]) : I in 1..X.len]).

"""

This prints

"""

[0,0,0] = 1
[0,0,1] = 1
[0,0,2] = 2
[0,0,3] = 3
[0,0,4] = 4
[0,1,0] = 1
[0,1,1] = 1
[0,1,2] = 2
[0,1,3] = 3
[0,1,4] = 4
[0,2,0] = 2
[0,2,1] = 2
[0,2,2] = 4
[0,2,3] = 6
[0,2,4] = 8
[0,3,0] = 3
[0,3,1] = 3
[0,3,2] = 6
[0,3,3] = 9
[0,3,4] = 12
[0,4,0] = 4
[0,4,1] = 4
[0,4,2] = 8
[0,4,3] = 12
[0,4,4] = 16
[1,0,0] = 1
[1,0,1] = 1
...

"""

Hope this solves your problem.

Best,

Hakan

[Nartoo Meon]

Hi, Hakan. Thank you!  (Why didn't cond/3 come to my head?)

" For a list of only 0s this will yield 1, which I assume is correct. "

Probably this doesn't affect the result. It's easy to add a constraint for no zeros at all positions once.

Background

Last 3 days I try to solve a weaker version of one "challenge" problem I introduce a year ago :

 

— Find an integer number that equal to power K of it's digits product. 

K must be > 1 and we need at least 2 digits for "interesting" solutions.

But, any number can be represented differently in different bases, so we may find solutions in other bases.

This may be the Only way to find a solution at all, because there are No solutions in base 10 up to number 10^305.
This is the result from two of my friends. First one wrote a program on Python to do some "intelligent" search even for some rational powers
and he proves that an infinite amount of solutions exist among all bases.
Second one proved that solutions exist only for powers K =< some irrational number between 6 and 7 . (formula with log)

My initial version restricted to power 2 instead of any possible K. So, here is some founded solutions for K=2:
13 in base 6, value: 9
22 in base 7, value: 16

523251 in base 7, value: 90_000

212323421441324231 in base 5, value: 2^28 * 3^8

(18 digits, the greatest they found, if I not lie)

Most solutions equal to a power of 2, especially for K>2 but only few of them have other prime factors like 3,5,7,11 and 17 (not sure).
An interesting solution exists in base 24, it's divisible by 23 and not so big. But, I forget it.

Bases 4 and 9 have the same behavior as 10. No solutions found (up to some number).

We see that solutions are rare.
Allowing an error ∓1 , we found only 3 solutions in base 10 up to 0.2 million.

So, to make solutions more common I need weaker conditions.
— Here is Why I need a product_except_0 constraint. 

My idea is to use polynomial representation, or even without a representation:

— Whitch polynomial P(X) with integer coefficients from 0..(Base–1) has a value P(Base) 

equal to product of it's coefficients, rised to a power K ?

But we often write polynomials like x^4 + x + 2 instead of  x^4 + 0*x^3 + 0*x^2 + x + 2
like number representation 10012, where 0's simply destroy the product.

Seems very natural.
All we need is that the coefficient against the highest power must be  > 0.

Possible problems

Currently, I don't know how to make a search across polynomials of different degrees at once, without calling the same predicate in a loop.
I can't add a new DVar into a list dynamically and use it later (at runtime), …or yes? 

Anyway I need an upper bound for degree, so I can use one (long) list of coefficients.

So, that's it. No more questions in my head now.

ср, 27 серп. 2025 р. о 08:49 Hakan Kjellerstrand <hak...@gmail.com> пише:

--
You received this message because you are subscribed to the Google Groups "Picat" group.
To unsubscribe from this group and stop receiving emails from it, send an email to picat-lang+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/picat-lang/02f7b51d-4ca3-479c-8d4a-f00470e5bdb3n%40googlegroups.com.

[Hakan Kjellerstrand]

Hi, Nartoo.

Thanks for the background on the product except 0. That's a fun problem!

Unless I'm missing something, here's a quite simple CP model for this:

"""

import cp. % Base 3..7: 1min33.5s
% import sat. % Base 3..7: 9min47.9s

main => go.

go ?=>
  nolog,
  member(Base,3..7),
  println(base=Base),
  member(Len,2..10), % number of digits
  member(K,2..5),  

  power_k(Base,K,Len, A,X) ,
  println([base=Base,len=Len,k=K,a=A,x_in_dec=X,x_in_base=X.to_radix_string(Base)]),
  fail, 
  nl.
go => true.

power_k(Base,K,Len, A,X) =>
  A = new_list(Len), % The list of digits
  A :: 0..Base-1,

  % No leading 0
  A[1] #!= 0,

  % convert to a number X (in base Base)
  to_num(A,Base,X),

  % Each digit ** K
  % Here is the product except 0 logic
  X #= prod([cond(T #= 0, 1,T)  : I in 1..Len, T #= A[I]**K]),

  solve($[ff,split],A ++ [X]).

% converts a number Num to/from a list of integer List given a base Base
to_num(List, Base, Num) =>
        Len = length(List),
        Num #= sum([List[I]*Base**(Len-I) : I in 1..Len]).

"""

The output is

"""

  base = 3
  [base = 3,len = 3,k = 4,a = [1,2,1],x_in_dec = 16,x_in_base = 121]
  [base = 3,len = 4,k = 5,a = [1,0,1,2],x_in_dec = 32,x_in_base = 1012]
  [base = 3,len = 7,k = 5,a = [1,1,0,1,2,2,1],x_in_dec = 1024,x_in_base = 1101221]
  [base = 3,len = 8,k = 4,a = [1,2,1,2,1,2,0,1],x_in_dec = 4096,x_in_base = 12121201]
  base = 4
  [base = 4,len = 2,k = 3,a = [2,0],x_in_dec = 8,x_in_base = 20]
  [base = 4,len = 3,k = 5,a = [2,0,0],x_in_dec = 32,x_in_base = 200]
  [base = 4,len = 4,k = 3,a = [3,1,2,0],x_in_dec = 216,x_in_base = 3120]
  base = 5
  [base = 5,len = 3,k = 4,a = [3,1,1],x_in_dec = 81,x_in_base = 311]
  [base = 5,len = 3,k = 5,a = [1,1,2],x_in_dec = 32,x_in_base = 112]
  [base = 5,len = 4,k = 2,a = [1,0,3,4],x_in_dec = 144,x_in_base = 1034]
  [base = 5,len = 9,k = 3,a = [2,1,1,3,0,2,4,2,1],x_in_dec = 884736,x_in_base = 211302421]
  base = 6
  [base = 6,len = 2,k = 2,a = [1,3],x_in_dec = 9,x_in_base = 13]
  [base = 6,len = 2,k = 3,a = [1,2],x_in_dec = 8,x_in_base = 12]
  [base = 6,len = 4,k = 2,a = [1,5,0,4],x_in_dec = 400,x_in_base = 1504]
  [base = 6,len = 4,k = 3,a = [2,2,1,2],x_in_dec = 512,x_in_base = 2212]
  [base = 6,len = 4,k = 4,a = [1,1,0,4],x_in_dec = 256,x_in_base = 1104]
  [base = 6,len = 6,k = 3,a = [3,2,5,0,0,0],x_in_dec = 27000,x_in_base = 325000]
  [base = 6,len = 6,k = 3,a = [4,1,1,4,1,2],x_in_dec = 32768,x_in_base = 411412]
  base = 7
  [base = 7,len = 2,k = 2,a = [2,2],x_in_dec = 16,x_in_base = 22]
  [base = 7,len = 6,k = 2,a = [5,2,3,2,5,1],x_in_dec = 90000,x_in_base = 523251]

"""

On my machine that took about 1min33s. On this specific test, the SAT solver is much slower (it took over 9 minutes), but it might be faster than the CP solver for larger instances.

I also did some other tests and found some more examples:

[base = 3,len = 12,k = 6,a = [1,1,1,0,2,2,1,2,1,0,0,1],x_in_dec = 262144,x_in_base = 111022121001]
[base = 3,len = 13,k = 5,a = [1,2,2,2,0,2,1,1,0,1,0,1,1],x_in_dec = 1048576,x_in_base = 1222021101011]
[base = 3,len = 14,k = 3,a = [1,0,2,2,1,1,1,2,2,0,2,0,2,2],x_in_dec = 2097152,x_in_base = 10221112202022]
[base = 3,len = 16,k = 5,a = [2,1,0,0,0,1,0,2,0,2,0,0,0,2,0,2],x_in_dec = 33554432,x_in_base = 2100010202000202]

[base = 5,len = 18,k = 2,a = [2,1,2,3,2,3,4,2,1,4,4,1,3,2,4,2,3,1],x_in_dec = 1761205026816,x_in_base = 212323421441324231]

Note that this CP model cannot handle numbers larger than 2**56 (base 10) since that's the limit of the constraint solvers.

But!

This is a change to play a little with the new bv module, introduced in the new version 3.9 which supports much larger domain variables.

Note that this use my bv_util.pi module found at https://hakank.org/picat/bv_utils.pi.

There are some caveats, however:

Unfortunately, right now I don't know  how to do the product_except_0 constraint for bv variables, so it will only give solutions without any 0s which is - of course - quite restricted.

It's also slower than the cp solver. 

But it can give solutions with very large values, so it might worth to explore it.

"""

import sat.
import bv.

import bv_utils. % https://hakank.org/picat/bv_utils.pi.

main => go.

go ?=>
  nolog,
  member(Base,3..10),
  println(base=Base),
  member(Len,2..20), % number of digits
  member(K,2..15),

  power_k_bv(Base,K,Len, A,X),
  println([base=Base,len=Len,k=K,a=A.bti,x_in_dec=X.bti,x_in_base=X.bti.to_radix_string(Base)]),
  fail,
  nl.
go => true.

power_k_bv(Base,K,Len, A,X) =>

  % A = make_bv_array(Len, 0, Base-1),
  A = make_bv_array(Len, 1, Base-1), % Skipping 0s

  % No leading 0
  bv_gt(A[1],0),

  % convert to a number X (in base Base)
  to_num_bv(A,Base,X),

  % Each digit ** K
  bv_eq(bv_prod([bv_pow(A[I],K) : I in 1..Len]),X), % Note: no 0 correction
 
  % product except 0 using cond/3 does not work
  % bv_eq(bv_prod([cond(T.bv_eq(0),1,T) : I in 1..Len, bv_pow(A[I],K).bv_eq(T) ]),X),

  solve(A ++ X).
"""

For the same setup as the cp model above, it gives these solutions (no numbers with 0s). This took about 2min47.3s (vs 1.33s using plain cp model above).

"""

base = 3
[base = 3,len = 3,k = 4,a = [1,2,1],x_in_dec = 16,x_in_base = 121]
base = 4
base = 5
[base = 5,len = 3,k = 4,a = [3,1,1],x_in_dec = 81,x_in_base = 311]
[base = 5,len = 3,k = 5,a = [1,1,2],x_in_dec = 32,x_in_base = 112]
base = 6
[base = 6,len = 2,k = 2,a = [1,3],x_in_dec = 9,x_in_base = 13]
[base = 6,len = 2,k = 3,a = [1,2],x_in_dec = 8,x_in_base = 12]
[base = 6,len = 4,k = 3,a = [2,2,1,2],x_in_dec = 512,x_in_base = 2212]
[base = 6,len = 6,k = 3,a = [4,1,1,4,1,2],x_in_dec = 32768,x_in_base = 411412]
base = 7
[base = 7,len = 2,k = 2,a = [2,2],x_in_dec = 16,x_in_base = 22]
[base = 7,len = 6,k = 2,a = [5,2,3,2,5,1],x_in_dec = 90000,x_in_base = 523251]
"""

Finding this takes about 25s:

[base = 5,len = 18,k = 2,a = [2,1,2,3,2,3,4,2,1,4,4,1,3,2,4,2,3,1],x_in_dec = 1761205026816,x_in_base = 212323421441324231]

As mentioned some days ago, I've done some more experiments with the bv module. see https://hakank.org/picat/#bv .

/Hakan

[David Silverman]

% X #= prod([cond(T #= 0, 1,T) : I in 1..Len, T #= A[I]**K]),

X #= prod([max(1,T) : I in 1..Len, T #= A[I]**K]),

 

I found that using max(1,T) instead of cond(T#=0,1,T) executes in about half the time when I tested it.