From mboxrd@z Thu Jan 1 00:00:00 1970 Path: news.gmane.io!.POSTED.blaine.gmane.org!not-for-mail From: Matt Wette Newsgroups: gmane.lisp.guile.user Subject: Re: Double-precision floating point arithmetics behaves differently in Guile and in C Date: Sun, 26 Sep 2021 12:11:03 -0700 Message-ID: References: Mime-Version: 1.0 Content-Type: text/plain; charset=utf-8; format=flowed Content-Transfer-Encoding: 8bit Injection-Info: ciao.gmane.io; posting-host="blaine.gmane.org:116.202.254.214"; logging-data="19571"; mail-complaints-to="usenet@ciao.gmane.io" User-Agent: Mozilla/5.0 (X11; Linux x86_64; rv:78.0) Gecko/20100101 Thunderbird/78.13.0 To: guile-user@gnu.org Original-X-From: guile-user-bounces+guile-user=m.gmane-mx.org@gnu.org Sun Sep 26 21:11:34 2021 Return-path: Envelope-to: guile-user@m.gmane-mx.org Original-Received: from lists.gnu.org ([209.51.188.17]) by ciao.gmane.io with esmtps (TLS1.2:ECDHE_RSA_AES_256_GCM_SHA384:256) (Exim 4.92) (envelope-from ) id 1mUZYj-0004vC-LP for guile-user@m.gmane-mx.org; Sun, 26 Sep 2021 21:11:33 +0200 Original-Received: from localhost ([::1]:60908 helo=lists1p.gnu.org) by lists.gnu.org with esmtp (Exim 4.90_1) (envelope-from ) id 1mUZYi-0007Hu-Kh for guile-user@m.gmane-mx.org; Sun, 26 Sep 2021 15:11:32 -0400 Original-Received: from eggs.gnu.org ([2001:470:142:3::10]:34194) by lists.gnu.org with esmtps (TLS1.2:ECDHE_RSA_AES_256_GCM_SHA384:256) (Exim 4.90_1) (envelope-from ) id 1mUZYN-0007FD-Di for guile-user@gnu.org; Sun, 26 Sep 2021 15:11:11 -0400 Original-Received: from mail-pl1-x62e.google.com ([2607:f8b0:4864:20::62e]:41496) by eggs.gnu.org with esmtps (TLS1.2:ECDHE_RSA_AES_128_GCM_SHA256:128) (Exim 4.90_1) (envelope-from ) id 1mUZYL-00029H-Kl for guile-user@gnu.org; Sun, 26 Sep 2021 15:11:11 -0400 Original-Received: by mail-pl1-x62e.google.com with SMTP id x8so7411180plv.8 for ; Sun, 26 Sep 2021 12:11:07 -0700 (PDT) DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=20210112; h=subject:to:references:from:message-id:date:user-agent:mime-version :in-reply-to:content-transfer-encoding:content-language; bh=op7R/WiJLtFodjapnrScqPFj0OIIyBfNisoXyyf3gCo=; b=WsT+UIl7WPbkJZXrz0hhEaFgL88/JVllQpGEsQvuP/dHr37USlTonGB61vF0GSG6h2 gm84sr3Ww4/B8c73x5ehaqplJzVA/ykGDM86OM7jhQj3/XdrcyokNCQVwkyD96rjDiYY q2Npqx1MErva/AryWDVna1/EA8fkgAPTMCT6zkdGMGbQaI46m3RY8iLkQIhHuU7IiHUj wG3Z2eFWTsLMgo3qhTnvHrwbAMbPf7CuhkvOpeFSzQpG2/3Oy075Qw3kHUBMWKCTWah4 7xeztoPECTvfJrhwHQYuP7+tn0hqSzVNwYoAu6ix0epHZVclUAU9tCY/uSB+Xw9POAXD mdjA== X-Google-DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=1e100.net; s=20210112; h=x-gm-message-state:subject:to:references:from:message-id:date :user-agent:mime-version:in-reply-to:content-transfer-encoding :content-language; bh=op7R/WiJLtFodjapnrScqPFj0OIIyBfNisoXyyf3gCo=; b=7UztoKrBJzWmE8v7HmRJu9uzWmFEW9qy+Q+iIEXG5ifh6vLzlS3yNYfRPuE/cXFvRs sUvDby80l3tVOTj2ccscZ1jyU2hH0lyYzVv+uwgfX18oNbxM1mSHmWqFSE4sTAsYbMKn M8f3JuhRFEtaLAnyZJd4u9PzydvEThYW8V++kCg+w0NqS+vmrSoUsrNnGM4RaN0hjeZJ eJyyZFuMOp1sMh3bj0hCcQY0CCALJITBdUzrcwU2Jt3rL4JMW3rwVdorvr8uLDoAXppa y0WLRlfylvM98j3Hd/XNRRgis6QKNt35tWuEo+HX9+qxFlj467GncwVKKaq+hJa8sJeR hV7g== X-Gm-Message-State: AOAM530UzesS7f0gY2utln7q+f+HStU6wHhkQ1AN31CAJMc+aHGx1Gjv EQ2wxNykVqTBn01CCHrY8dTzD/fIn/c= X-Google-Smtp-Source: ABdhPJxxW2uU5xiCDBUGlKUFbxq+x4T1ISGcVxKQKubRsCt77Mnz3MRTrcHVuXKUk6T7PKY63xarsQ== X-Received: by 2002:a17:903:49:b0:13a:752e:244c with SMTP id l9-20020a170903004900b0013a752e244cmr18969660pla.81.1632683465643; Sun, 26 Sep 2021 12:11:05 -0700 (PDT) Original-Received: from [192.168.2.163] (64-52-176-132.championbroadband.com. [64.52.176.132]) by smtp.gmail.com with ESMTPSA id v6sm14571415pfv.83.2021.09.26.12.11.04 for (version=TLS1_3 cipher=TLS_AES_128_GCM_SHA256 bits=128/128); Sun, 26 Sep 2021 12:11:05 -0700 (PDT) In-Reply-To: Content-Language: en-US Received-SPF: pass client-ip=2607:f8b0:4864:20::62e; envelope-from=matt.wette@gmail.com; helo=mail-pl1-x62e.google.com X-Spam_score_int: -5 X-Spam_score: -0.6 X-Spam_bar: / X-Spam_report: (-0.6 / 5.0 requ) BAYES_00=-1.9, DKIM_SIGNED=0.1, DKIM_VALID=-0.1, DKIM_VALID_AU=-0.1, DKIM_VALID_EF=-0.1, FREEMAIL_FROM=0.001, NICE_REPLY_A=-0.478, RCVD_IN_DNSWL_NONE=-0.0001, SPF_HELO_NONE=0.001, SPF_PASS=-0.001, URI_DOTEDU=1.999 autolearn=no autolearn_force=no X-Spam_action: no action X-BeenThere: guile-user@gnu.org X-Mailman-Version: 2.1.23 Precedence: list List-Id: General Guile related discussions List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , Errors-To: guile-user-bounces+guile-user=m.gmane-mx.org@gnu.org Original-Sender: "guile-user" Xref: news.gmane.io gmane.lisp.guile.user:17784 Archived-At: Your post is a little vague.  Are you comparing singular values vectors of A vs eigenvalues, eigenvectors of A'*A ?  Also, SVD is iterative and different algorithms may yield different results (and some of the SVD algorithms out there are not great).  Are you using the exact same algorithm but noticing that (op x y) in Guile behaves differently than (x op y) in C? On 9/26/21 1:45 AM, Panicz Maciej Godek wrote: > Hi, > > I recently ran into a very weird problem. I am using code for computing > Singular Value Decomposition that's been ported do Scheme by Gerald Sussman > as a part of the scmutils package (more specifically, I use some bits of > code extracted from guile-scmutils > ) > > The code is too long to paste here, but I have uploaded a copy to my github > > . > > I have found an ellipsoid-fitting algorithm that leans on this SVD code: > when I run it on a randomly generated set of points lying on an ellipsoid, > the ellipsoid parameters are reconstructed correctly. > > The trouble begins when I try to port that code to C. Initially, I have > used a C code for computing SVD that I've found elsewhere (which was based > on the same Fortran code from LINPACK), but I eventually wrote code to > convert Scheme code to C to make sure that both algorithms are identical > (and at this point I am quite confident that they are). > > The problem is, that the results of floating-point operations behave > slightly differently in C than in Guile, and these small numerical > differences accumulate, and eventually start making huge differences -- at > least that's something that I've noticed. Moreover, I can quite confidently > say that the results that I get in C are incorrect -- that the *eigenvalues* > (which translate to the radii of the ellipsoid) are very different, the > calculated center of the ellipsoid is slightly different, and that the > *eigenvectors* (which can be interpreted as a rotation of the ellipsoid) > are very similar. > > Either way, the ellipsoid isn't reconstructed correctly, because the points > used for reconstruction do not lie near that ellipsoid (unlike in the case > of the C code). > > I am really troubled with this behavior -- according to the documentation, > and to the bits of source code that I managed to analyze, Guile uses > exactly the same representation of floating point numbers as the one I use > in C -- namely, the double type. (I only use single-precision floats to > express the input points, but I ran some tests using small integers on > inputs, to make sure, that they are exactly representable -- and the > results were still very different) > > Is there a way to make sure that the floating point operations in C behave > identically as those in Guile?